Index: /palm/trunk/SOURCE/Makefile
===================================================================
--- /palm/trunk/SOURCE/Makefile	(revision 3434)
+++ /palm/trunk/SOURCE/Makefile	(revision 3435)
@@ -25,4 +25,8 @@
 # -----------------
 # $Id$
+# - Add surface_mod to data_output_mask
+# - Add chemistry_model_mod and surface_mod to init_masks
+# 
+# 3421 2018-10-24 18:39:32Z gronemeier
 # Add netcdf_data_input_mod to netcdf_interface_mod
 # bugfix: add dependencies to chemistry_model_mod
@@ -879,5 +883,6 @@
 	mod_particle_attributes.o \
 	modules.o \
-	netcdf_interface_mod.o
+	netcdf_interface_mod.o \
+	surface_mod.o
 data_output_profiles.o: \
 	cpulog_mod.o \
@@ -1087,7 +1092,9 @@
 init_masks.o: \
 	bulk_cloud_model_mod.o \
-	mod_kinds.o \
-	modules.o \
-	netcdf_interface_mod.o
+	chemistry_model_mod.o \
+	mod_kinds.o \
+	modules.o \
+	netcdf_interface_mod.o \
+	radiation_model_mod.o
 init_pegrid.o: \
 	mod_kinds.o \
Index: /palm/trunk/SOURCE/chemistry_model_mod.f90
===================================================================
--- /palm/trunk/SOURCE/chemistry_model_mod.f90	(revision 3434)
+++ /palm/trunk/SOURCE/chemistry_model_mod.f90	(revision 3435)
@@ -27,4 +27,7 @@
 ! -----------------
 ! $Id$
+! Add terrain-following masked output
+! 
+! 3373 2018-10-18 15:25:56Z kanani
 ! Remove MPI_Abort, replace by message
 !
@@ -1149,7 +1152,9 @@
        USE kinds
        USE pegrid,             ONLY: myid, threads_per_task
-
+       USE surface_mod,        ONLY: get_topography_top_index_ji
 
        IMPLICIT NONE
+
+       CHARACTER(LEN=5) ::  grid        !< flag to distinquish between staggered grids
 
        CHARACTER (LEN=*)::  variable    !< 
@@ -1166,5 +1171,8 @@
        INTEGER(iwp) ::  j               !< grid index along y-direction
        INTEGER(iwp) ::  k               !< grid index along z-direction
+       INTEGER(iwp) ::  topo_top_ind    !< k index of highest horizontal surface
+
        found = .TRUE.
+       grid  = 's'
 
        spec_name = TRIM( variable(4:) )
@@ -1178,22 +1186,81 @@
 !                                                           TRIM(chem_species(lsp)%name)       
              IF (av == 0) THEN
-                DO  i = 1, mask_size_l(mid,1)
-                   DO  j = 1, mask_size_l(mid,2) 
-                      DO  k = 1, mask_size(mid,3) 
-                          local_pf(i,j,k) = chem_species(lsp)%conc(mask_k(mid,k),                                    &
-                                              mask_j(mid,j), mask_i(mid,i))
+                IF ( .NOT. mask_surface(mid) )  THEN
+
+                   DO  i = 1, mask_size_l(mid,1)
+                      DO  j = 1, mask_size_l(mid,2) 
+                         DO  k = 1, mask_size(mid,3) 
+                             local_pf(i,j,k) = chem_species(lsp)%conc(  &
+                                                  mask_k(mid,k),        &
+                                                  mask_j(mid,j),        &
+                                                  mask_i(mid,i)      )
+                         ENDDO
                       ENDDO
                    ENDDO
-                ENDDO
-           
-             ELSE
-                DO  i = 1, mask_size_l(mid,1)
-                   DO  j = 1, mask_size_l(mid,2)
-                      DO  k =  1, mask_size_l(mid,3)
-                          local_pf(i,j,k) = chem_species(lsp)%conc_av(mask_k(mid,k),                &
-                                             mask_j(mid,j), mask_i(mid,i))
+
+                ELSE
+!              
+!--                Terrain-following masked output
+                   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
+                         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)      )
+                         ENDDO
                       ENDDO
                    ENDDO
-                ENDDO
+
+                ENDIF
+             ELSE
+                IF ( .NOT. mask_surface(mid) )  THEN
+
+                   DO  i = 1, mask_size_l(mid,1)
+                      DO  j = 1, mask_size_l(mid,2)
+                         DO  k =  1, mask_size_l(mid,3)
+                             local_pf(i,j,k) = chem_species(lsp)%conc_av(  &
+                                                  mask_k(mid,k),           &
+                                                  mask_j(mid,j),           &
+                                                  mask_i(mid,i)         )
+                         ENDDO
+                      ENDDO
+                   ENDDO
+
+                ELSE
+!              
+!--                Terrain-following masked output
+                   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
+                         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)         )
+                         ENDDO
+                      ENDDO
+                   ENDDO
+
+                ENDIF
+
+
              ENDIF
              found = .FALSE.
Index: /palm/trunk/SOURCE/data_output_mask.f90
===================================================================
--- /palm/trunk/SOURCE/data_output_mask.f90	(revision 3434)
+++ /palm/trunk/SOURCE/data_output_mask.f90	(revision 3435)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Add terrain-following output
+! 
+! 3421 2018-10-24 18:39:32Z gronemeier
 ! Renamed output variables
 ! 
@@ -155,4 +158,5 @@
         ONLY:  air_chemistry, domask, domask_no, domask_time_count, mask_i,    &
                mask_j, mask_k, mask_size, mask_size_l, mask_start_l,           &
+               mask_surface,                                                   &
                max_masks, message_string, mid, nz_do3d, simulated_time
     USE cpulog,                                                                &
@@ -160,5 +164,5 @@
 
     USE indices,                                                               &
-        ONLY:  nbgp, nxl, nxr, nyn, nys, nzb
+        ONLY:  nbgp, nxl, nxr, nyn, nys, nzb, nzt
 
     USE kinds
@@ -182,23 +186,30 @@
         ONLY:  radiation, radiation_data_output_mask
 
+    USE surface_mod,                                                           &
+        ONLY :  surf_def_h, surf_lsm_h, surf_usm_h, get_topography_top_index_ji
+
     IMPLICIT NONE
 
-    INTEGER(iwp) ::  av       !< 
-    INTEGER(iwp) ::  ngp      !< 
-    INTEGER(iwp) ::  i        !< 
-    INTEGER(iwp) ::  ivar     !<
-    INTEGER(iwp) ::  j        !< 
-    INTEGER(iwp) ::  k        !< 
-    INTEGER(iwp) ::  n        !< 
+    CHARACTER(LEN=5) ::  grid !< flag to distinquish between staggered grids
+
+    INTEGER(iwp) ::  av                      !< 
+    INTEGER(iwp) ::  ngp                     !< 
+    INTEGER(iwp) ::  i                       !< 
+    INTEGER(iwp) ::  ivar                    !<
+    INTEGER(iwp) ::  j                       !< 
+    INTEGER(iwp) ::  k                       !< 
+    INTEGER(iwp) ::  kk                      !< 
+    INTEGER(iwp) ::  n                       !< 
     INTEGER(iwp) ::  netcdf_data_format_save !<
-    INTEGER(iwp) ::  sender   !< 
-    INTEGER(iwp) ::  ind(6)   !< 
-    
-    LOGICAL ::  found         !< 
-    LOGICAL ::  resorted      !< 
-    
-    REAL(wp) ::  mean_r       !< 
-    REAL(wp) ::  s_r2         !< 
-    REAL(wp) ::  s_r3         !< 
+    INTEGER(iwp) ::  sender                  !< 
+    INTEGER(iwp) ::  topo_top_ind            !< k index of highest horizontal surface
+    INTEGER(iwp) ::  ind(6)                  !< 
+
+    LOGICAL ::  found     !< 
+    LOGICAL ::  resorted  !< 
+
+    REAL(wp) ::  mean_r   !< 
+    REAL(wp) ::  s_r2     !< 
+    REAL(wp) ::  s_r3     !< 
     
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  local_pf    !<
@@ -261,4 +272,7 @@
        ENDIF
 !
+!--    Set default grid for terrain-following output
+       grid = 's'
+!
 !--    Set flag to steer output of radiation, land-surface, or user-defined
 !--    quantities
@@ -308,12 +322,32 @@
                 tend = prt_count
                 CALL exchange_horiz( tend, nbgp )
-                DO  i = 1, mask_size_l(mid,1)
-                   DO  j = 1, mask_size_l(mid,2)
-                      DO  k = 1, mask_size_l(mid,3)
-                         local_pf(i,j,k) =  tend(mask_k(mid,k), &
-                                   mask_j(mid,j),mask_i(mid,i))
-                      ENDDO
-                   ENDDO
-                ENDDO
+                IF ( .NOT. mask_surface(mid) )  THEN
+                   DO  i = 1, mask_size_l(mid,1)
+                      DO  j = 1, mask_size_l(mid,2)
+                         DO  k = 1, mask_size_l(mid,3)
+                            local_pf(i,j,k) =  tend(mask_k(mid,k), &
+                                      mask_j(mid,j),mask_i(mid,i))
+                         ENDDO
+                      ENDDO
+                   ENDDO
+                ELSE
+!
+!--                Terrain-following masked output
+                   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 )
+                         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))
+                         ENDDO
+                      ENDDO
+                   ENDDO
+                ENDIF
                 resorted = .TRUE.
              ELSE
@@ -354,12 +388,32 @@
                    tend = 0.0_wp
                 ENDIF
-                DO  i = 1, mask_size_l(mid,1)
-                   DO  j = 1, mask_size_l(mid,2)
-                      DO  k = 1, mask_size_l(mid,3)
-                         local_pf(i,j,k) =  tend(mask_k(mid,k), &
-                                   mask_j(mid,j),mask_i(mid,i))
-                      ENDDO
-                   ENDDO
-                ENDDO
+                IF ( .NOT. mask_surface(mid) )  THEN
+                   DO  i = 1, mask_size_l(mid,1)
+                      DO  j = 1, mask_size_l(mid,2)
+                         DO  k = 1, mask_size_l(mid,3)
+                            local_pf(i,j,k) =  tend(mask_k(mid,k), &
+                                      mask_j(mid,j),mask_i(mid,i))
+                         ENDDO
+                      ENDDO
+                   ENDDO
+                ELSE
+!
+!--                Terrain-following masked output
+                   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 )
+                         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))
+                         ENDDO
+                      ENDDO
+                   ENDDO
+                ENDIF
                 resorted = .TRUE.
              ELSE
@@ -373,14 +427,36 @@
                    to_be_resorted => pt
                 ELSE
-                   DO  i = 1, mask_size_l(mid,1)
-                      DO  j = 1, mask_size_l(mid,2)
-                         DO  k = 1, mask_size_l(mid,3)
-                            local_pf(i,j,k) =  &
-                                 pt(mask_k(mid,k),mask_j(mid,j),mask_i(mid,i)) &
-                                 + lv_d_cp * d_exner(mask_k(mid,k)) *          &
-                                   ql(mask_k(mid,k),mask_j(mid,j),mask_i(mid,i))
-                         ENDDO
-                      ENDDO
-                   ENDDO
+                   IF ( .NOT. mask_surface(mid) )  THEN
+                      DO  i = 1, mask_size_l(mid,1)
+                         DO  j = 1, mask_size_l(mid,2)
+                            DO  k = 1, mask_size_l(mid,3)
+                               local_pf(i,j,k) =  &
+                                  pt(mask_k(mid,k),mask_j(mid,j),mask_i(mid,i)) &
+                                  + lv_d_cp * d_exner(mask_k(mid,k)) *          &
+                                    ql(mask_k(mid,k),mask_j(mid,j),mask_i(mid,i))
+                            ENDDO
+                         ENDDO
+                      ENDDO
+                   ELSE
+!
+!--                   Terrain-following masked output
+                      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 )
+                            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))
+                            ENDDO
+                         ENDDO
+                      ENDDO
+                   ENDIF                   
                    resorted = .TRUE.
                 ENDIF
@@ -447,12 +523,32 @@
                    tend = 0.0_wp
                 ENDIF
-                DO  i = 1, mask_size_l(mid,1)
-                   DO  j = 1, mask_size_l(mid,2)
-                      DO  k = 1, mask_size_l(mid,3)
-                         local_pf(i,j,k) =  tend(mask_k(mid,k), &
-                                   mask_j(mid,j),mask_i(mid,i))
-                      ENDDO
-                   ENDDO
-                ENDDO
+                IF ( .NOT. mask_surface(mid) )  THEN
+                   DO  i = 1, mask_size_l(mid,1)
+                      DO  j = 1, mask_size_l(mid,2)
+                         DO  k = 1, mask_size_l(mid,3)
+                            local_pf(i,j,k) =  tend(mask_k(mid,k), &
+                                      mask_j(mid,j),mask_i(mid,i))
+                         ENDDO
+                      ENDDO
+                   ENDDO
+                ELSE
+!
+!--                Terrain-following masked output
+                   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 )
+                         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))
+                         ENDDO
+                      ENDDO
+                   ENDDO
+                ENDIF
                 resorted = .TRUE.
              ELSE
@@ -463,13 +559,34 @@
           CASE ( 'qv' )
              IF ( av == 0 )  THEN
-                DO  i = 1, mask_size_l(mid,1)
-                   DO  j = 1, mask_size_l(mid,2)
-                      DO  k = 1, mask_size_l(mid,3)
-                         local_pf(i,j,k) =  &
-                              q(mask_k(mid,k),mask_j(mid,j),mask_i(mid,i)) -  &
-                              ql(mask_k(mid,k),mask_j(mid,j),mask_i(mid,i))
-                      ENDDO
-                   ENDDO
-                ENDDO
+                IF ( .NOT. mask_surface(mid) )  THEN
+                   DO  i = 1, mask_size_l(mid,1)
+                      DO  j = 1, mask_size_l(mid,2)
+                         DO  k = 1, mask_size_l(mid,3)
+                            local_pf(i,j,k) =  &
+                                 q(mask_k(mid,k),mask_j(mid,j),mask_i(mid,i)) -  &
+                                 ql(mask_k(mid,k),mask_j(mid,j),mask_i(mid,i))
+                         ENDDO
+                      ENDDO
+                   ENDDO
+                ELSE
+!
+!--                Terrain-following masked output
+                   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 )
+                         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))
+                         ENDDO
+                      ENDDO
+                   ENDDO
+                ENDIF
                 resorted = .TRUE.
              ELSE
@@ -527,4 +644,5 @@
 
           CASE ( 'w' )
+             grid = 'w'
              IF ( av == 0 )  THEN
                 to_be_resorted => w
@@ -566,12 +684,39 @@
 !--    Resort the array to be output, if not done above
        IF ( .NOT. resorted )  THEN
-          DO  i = 1, mask_size_l(mid,1)
-             DO  j = 1, mask_size_l(mid,2)
-                DO  k = 1, mask_size_l(mid,3)
-                   local_pf(i,j,k) =  to_be_resorted(mask_k(mid,k), &
-                                      mask_j(mid,j),mask_i(mid,i))
+          IF ( .NOT. mask_surface(mid) )  THEN
+!
+!--          Default masked output
+             DO  i = 1, mask_size_l(mid,1)
+                DO  j = 1, mask_size_l(mid,2)
+                   DO  k = 1, mask_size_l(mid,3)
+                      local_pf(i,j,k) =  to_be_resorted(mask_k(mid,k), &
+                                         mask_j(mid,j),mask_i(mid,i))
+                   ENDDO
                 ENDDO
              ENDDO
-          ENDDO
+
+          ELSE
+!
+!--          Terrain-following masked output
+             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
+                   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)                     )
+                   ENDDO
+                ENDDO
+             ENDDO
+
+          ENDIF
        ENDIF
 
@@ -710,3 +855,4 @@
 #endif
 
+
  END SUBROUTINE data_output_mask
Index: /palm/trunk/SOURCE/init_masks.f90
===================================================================
--- /palm/trunk/SOURCE/init_masks.f90	(revision 3434)
+++ /palm/trunk/SOURCE/init_masks.f90	(revision 3435)
@@ -25,4 +25,8 @@
 ! -----------------
 ! $Id$
+! Add checks for chemistry and radiation model
+! Set indices for terrain-following masked output
+! 
+! 3421 2018-10-24 18:39:32Z gronemeier
 ! Renamed output variables
 ! 
@@ -142,16 +146,24 @@
         ONLY:  zu, zw
 
+    USE bulk_cloud_model_mod,                                                  &
+        ONLY: bulk_cloud_model, microphysics_morrison, microphysics_seifert
+
+    USE chemistry_model_mod,                                                   &
+        ONLY:  chem_check_data_output
+
     USE control_parameters,                                                    &
-        ONLY:  constant_diffusion, cloud_droplets,                             &
+        ONLY:  air_chemistry,                                                  &
+               constant_diffusion, cloud_droplets,                             &
                data_output_masks, data_output_masks_user,                      &
                doav, doav_n, domask, domask_no, dz, dz_stretch_level_start,    &
                humidity, mask, masks, mask_scale, mask_i,                      &
                mask_i_global, mask_j, mask_j_global, mask_k, mask_k_global,    &
-               mask_loop, mask_size, mask_size_l, mask_start_l, mask_x,        &
+               mask_k_over_surface,                                            &
+               mask_loop, mask_size, mask_size_l, mask_start_l,                &
+               mask_surface, mask_x,                                           &
                mask_x_loop, mask_xyz_dimension, mask_y, mask_y_loop, mask_z,   &
                mask_z_loop, max_masks,  message_string, mid,                   &
                passive_scalar, ocean_mode, varnamelength
 
-
     USE grid_variables,                                                        &
         ONLY:  dx, dy
@@ -162,7 +174,4 @@
     USE kinds
 
-    USE bulk_cloud_model_mod,                                                  &
-        ONLY: bulk_cloud_model, microphysics_morrison, microphysics_seifert
-
     USE netcdf_interface,                                                      &
         ONLY:  domask_unit, netcdf_data_format
@@ -172,4 +181,7 @@
 
     USE pegrid
+
+    USE radiation_model_mod,                                                   &
+        ONLY:  radiation, radiation_check_data_output
 
     IMPLICIT NONE
@@ -181,4 +193,5 @@
     CHARACTER (LEN=varnamelength), DIMENSION(max_masks,100) ::  do_mask_user !<
 
+    INTEGER(iwp) ::  count        !< counting masking indices along a dimension
     INTEGER(iwp) ::  i            !< 
     INTEGER(iwp) ::  ilen         !<
@@ -187,4 +200,5 @@
     INTEGER(iwp) ::  j            !<
     INTEGER(iwp) ::  k            !<
+    INTEGER(iwp) ::  m            !< mask index
     INTEGER(iwp) ::  n            !<
     INTEGER(iwp) ::  sender       !<
@@ -232,5 +246,10 @@
        mask      (mid,2,:) = mask_y(mid,:)
        mask      (mid,3,:) = mask_z(mid,:) 
-       
+!
+!--    Flag a mask as terrain following
+       IF ( mask_k_over_surface(mid,1) /= -1_iwp )  THEN
+          mask_surface(mid) = .TRUE.
+       ENDIF
+
        IF ( mask_x_loop(mid,1) == -1.0_wp  .AND.  mask_x_loop(mid,2) == -1.0_wp&
             .AND.  mask_x_loop(mid,3) == -1.0_wp )  THEN
@@ -456,5 +475,15 @@
 
              CASE DEFAULT
+             
                 CALL user_check_data_output( var, unit )
+
+                IF ( unit == 'illegal'  .AND.  air_chemistry                   &
+                     .AND.  (var(1:3) == 'kc_' .OR. var(1:3) == 'em_') )  THEN
+                   CALL chem_check_data_output( var, unit, 0, 0, 0 )
+                ENDIF
+
+                IF ( unit == 'illegal' )  THEN
+                   CALL radiation_check_data_output( var, unit, 0, 0, 0 )
+                ENDIF
 
                 IF ( unit == 'illegal' )  THEN
@@ -505,5 +534,27 @@
        CALL set_mask_locations( 1, dx, 'dx', nx, 'nx', nxl, nxr )
        CALL set_mask_locations( 2, dy, 'dy', ny, 'ny', nys, nyn )
-       CALL set_mask_locations( 3, dz(1), 'dz', nz, 'nz', nzb, nzt )
+       IF ( .NOT. mask_surface(mid) )  THEN
+          CALL set_mask_locations( 3, dz(1), 'dz', nz, 'nz', nzb, nzt )
+       ELSE
+!
+!--       Set vertical mask locations and size in case of terrain-following
+!--       output
+          count = 0
+          DO  WHILE ( mask_k_over_surface(mid, count+1) >= 0 )
+             m = mask_k_over_surface(mid, count+1)
+             IF ( m > nz+1 )  THEN
+                WRITE ( message_string, '(I3,A,I3,A,I1,3A,I3)' )               &
+                     m,' in mask ',mid,' along dimension ', 3,                 &
+                     ' exceeds (nz+1) = ',nz+1
+                CALL message( 'init_masks', 'PA0331', 1, 2, 0, 6, 0 )
+             ENDIF
+             count = count + 1
+             mask_k(mid,count) = mask_k_over_surface(mid, count)
+             IF ( count == mask_xyz_dimension )  EXIT
+          ENDDO
+          mask_start_l(mid,3) = 1
+          mask_size(mid,3)    = count
+          mask_size_l(mid,3)  = count
+       ENDIF
 !
 !--    Set global masks along all three dimensions (required by
Index: /palm/trunk/SOURCE/modules.f90
===================================================================
--- /palm/trunk/SOURCE/modules.f90	(revision 3434)
+++ /palm/trunk/SOURCE/modules.f90	(revision 3435)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! +mask_k_over_surface, mask_surface
+! 
+! 3422 2018-10-24 19:01:57Z gronemeier
 ! bugfix: increase number of blanks in output string
 ! 
@@ -1283,4 +1286,6 @@
     INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  mask_j_global  !< global grid index of masked output point on y-dimension
     INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  mask_k_global  !< global grid index of masked output point on z-dimension
+
+    INTEGER(iwp), DIMENSION(max_masks,mask_xyz_dimension) ::  mask_k_over_surface = -1  !< namelist parameter, k index of height over surface
 
     LOGICAL ::  agent_time_unlimited = .FALSE.                   !< namelist parameter
@@ -1393,4 +1398,6 @@
     LOGICAL ::  data_output_xz(0:1) = .FALSE.                !< output of xz cross-section data?
     LOGICAL ::  data_output_yz(0:1) = .FALSE.                !< output of yz cross-section data?
+
+    LOGICAL, DIMENSION(max_masks) ::  mask_surface = .FALSE.      !< flag for surface-following masked output
 
     REAL(wp) ::  advected_distance_x = 0.0_wp                  !< advected distance of model domain along x 
Index: /palm/trunk/SOURCE/netcdf_interface_mod.f90
===================================================================
--- /palm/trunk/SOURCE/netcdf_interface_mod.f90	(revision 3434)
+++ /palm/trunk/SOURCE/netcdf_interface_mod.f90	(revision 3435)
@@ -25,4 +25,8 @@
 ! -----------------
 ! $Id$
+! Bugfix: corrected order of calls to define_netcdf_grid for masked output
+! Add vertical dimensions to masked output in case of terrain-following output
+! 
+! 3421 2018-10-24 18:39:32Z gronemeier
 ! Bugfix: move ocean output variables to ocean_mod
 ! Renamed output variables
@@ -589,5 +593,6 @@
                do3d_time_count, domask_time_count, end_time, land_surface,     &
                mask_size_l, mask_i, mask_i_global, mask_j, mask_j_global,      &
-               mask_k_global, message_string, mid, ntdim_2d_xy, ntdim_2d_xz,   &
+               mask_k_global, mask_surface,                                    &
+               message_string, mid, ntdim_2d_xy, ntdim_2d_xz,                  &
                ntdim_2d_yz, ntdim_3d, nz_do3d, ocean_mode, plant_canopy,       &
                run_description_header, section, simulated_time,                &
@@ -873,21 +878,46 @@
 !
 !--       Define spatial dimensions and coordinates:
-!--       Define vertical coordinate grid (zu grid)
-          CALL netcdf_create_dim( id_set_mask(mid,av), 'zu_3d',                &
-                                  mask_size(mid,3), id_dim_zu_mask(mid,av),    &
-                                  470 )
-          CALL netcdf_create_var( id_set_mask(mid,av),                         &
-                                  (/ id_dim_zu_mask(mid,av) /), 'zu_3d',       &
-                                  NF90_DOUBLE, id_var_zu_mask(mid,av),         &
-                                  'meters', '', 471, 472, 000 )
-!
-!--       Define vertical coordinate grid (zw grid)
-          CALL netcdf_create_dim( id_set_mask(mid,av), 'zw_3d',                &
-                                  mask_size(mid,3), id_dim_zw_mask(mid,av),    &
-                                  473 )
-          CALL netcdf_create_var( id_set_mask(mid,av),                         &
-                                  (/ id_dim_zw_mask(mid,av) /), 'zw_3d',       &
-                                  NF90_DOUBLE, id_var_zw_mask(mid,av),         &
-                                 'meters', '', 474, 475, 000 )
+          IF ( mask_surface(mid) )  THEN
+!
+!--          In case of terrain-following output, the vertical dimensions are
+!--          indices, not meters
+             CALL netcdf_create_dim( id_set_mask(mid,av), 'ku_above_surf',     &
+                                     mask_size(mid,3), id_dim_zu_mask(mid,av), &
+                                     470 )
+             CALL netcdf_create_var( id_set_mask(mid,av),                      &
+                                     (/ id_dim_zu_mask(mid,av) /),             &
+                                     'ku_above_surf',                          &
+                                     NF90_DOUBLE, id_var_zu_mask(mid,av),      &
+                                     '1', 'grid point above terrain',          &
+                                     471, 472, 000 )
+             CALL netcdf_create_dim( id_set_mask(mid,av), 'kw_above_surf',     &
+                                     mask_size(mid,3), id_dim_zw_mask(mid,av), &
+                                     473 )
+             CALL netcdf_create_var( id_set_mask(mid,av),                      &
+                                     (/ id_dim_zw_mask(mid,av) /),             &
+                                     'kw_above_surf',                          &
+                                     NF90_DOUBLE, id_var_zw_mask(mid,av),      &
+                                    '1', 'grid point above terrain',           &
+                                    474, 475, 000 )
+          ELSE
+!
+!--          Define vertical coordinate grid (zu grid)
+             CALL netcdf_create_dim( id_set_mask(mid,av), 'zu_3d',             &
+                                     mask_size(mid,3), id_dim_zu_mask(mid,av), &
+                                     470 )
+             CALL netcdf_create_var( id_set_mask(mid,av),                      &
+                                     (/ id_dim_zu_mask(mid,av) /), 'zu_3d',    &
+                                     NF90_DOUBLE, id_var_zu_mask(mid,av),      &
+                                     'meters', '', 471, 472, 000 )
+!
+!--          Define vertical coordinate grid (zw grid)
+             CALL netcdf_create_dim( id_set_mask(mid,av), 'zw_3d',             &
+                                     mask_size(mid,3), id_dim_zw_mask(mid,av), &
+                                     473 )
+             CALL netcdf_create_var( id_set_mask(mid,av),                      &
+                                     (/ id_dim_zw_mask(mid,av) /), 'zw_3d',    &
+                                     NF90_DOUBLE, id_var_zw_mask(mid,av),      &
+                                    'meters', '', 474, 475, 000 )
+          ENDIF
 !
 !--       Define x-axis (for scalar position)
@@ -1068,4 +1098,7 @@
 !
 !--                Check for quantities defined in other modules
+                   CALL tcm_define_netcdf_grid( domask( mid,av,i), found,      &
+                                                        grid_x, grid_y, grid_z )
+
                    IF ( .NOT. found  .AND.  air_chemistry )  THEN
                       CALL chem_define_netcdf_grid( domask(mid,av,i), found,   &
@@ -1105,8 +1138,4 @@
                                                          grid_z )
                    ENDIF
-
-                   CALL tcm_define_netcdf_grid( domask( mid,av,i), found,      &
-                                                        grid_x, grid_y, grid_z )
-
 !
 !--                Now check for user-defined quantities
@@ -1343,17 +1372,37 @@
           ALLOCATE( netcdf_data(mask_size(mid,3)) )
 
-          netcdf_data = zu( mask_k_global(mid,:mask_size(mid,3)) )
-
-          nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), id_var_zu_mask(mid,av), &
-                                  netcdf_data, start = (/ 1 /), &
-                                  count = (/ mask_size(mid,3) /) )
-          CALL netcdf_handle_error( 'netcdf_define_header', 503 )
-
-          netcdf_data = zw( mask_k_global(mid,:mask_size(mid,3)) )
-
-          nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), id_var_zw_mask(mid,av), &
-                                  netcdf_data, start = (/ 1 /), &
-                                  count = (/ mask_size(mid,3) /) )
-          CALL netcdf_handle_error( 'netcdf_define_header', 504 )
+          IF ( mask_surface(mid) )  THEN
+
+             netcdf_data = mask_k_global(mid,:mask_size(mid,3))
+
+             nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), id_var_zu_mask(mid,av), &
+                                     netcdf_data, start = (/ 1 /), &
+                                     count = (/ mask_size(mid,3) /) )
+             CALL netcdf_handle_error( 'netcdf_define_header', 503 )
+
+             netcdf_data = mask_k_global(mid,:mask_size(mid,3))
+
+             nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), id_var_zw_mask(mid,av), &
+                                     netcdf_data, start = (/ 1 /), &
+                                     count = (/ mask_size(mid,3) /) )
+             CALL netcdf_handle_error( 'netcdf_define_header', 504 )
+
+          ELSE
+
+             netcdf_data = zu( mask_k_global(mid,:mask_size(mid,3)) )
+
+             nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), id_var_zu_mask(mid,av), &
+                                     netcdf_data, start = (/ 1 /), &
+                                     count = (/ mask_size(mid,3) /) )
+             CALL netcdf_handle_error( 'netcdf_define_header', 503 )
+
+             netcdf_data = zw( mask_k_global(mid,:mask_size(mid,3)) )
+
+             nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), id_var_zw_mask(mid,av), &
+                                     netcdf_data, start = (/ 1 /), &
+                                     count = (/ mask_size(mid,3) /) )
+             CALL netcdf_handle_error( 'netcdf_define_header', 504 )
+
+          ENDIF
 
           DEALLOCATE( netcdf_data )
Index: /palm/trunk/SOURCE/parin.f90
===================================================================
--- /palm/trunk/SOURCE/parin.f90	(revision 3434)
+++ /palm/trunk/SOURCE/parin.f90	(revision 3435)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Add mask_k_over_surface
+! 
+! 3421 2018-10-24 18:39:32Z gronemeier
 ! Added module for data output at surfaces
 ! 
@@ -698,5 +701,6 @@
              dt_domask, dt_dopr, dt_dopr_listing, dt_dots, dt_do2d_xy,         &
              dt_do2d_xz, dt_do2d_yz, dt_do3d, dt_max, dt_restart,              &
-             dt_run_control,end_time, force_print_header, mask_scale_x,        &
+             dt_run_control,end_time, force_print_header, mask_k_over_surface, &
+             mask_scale_x,        &
              mask_scale_y, mask_scale_z, mask_x, mask_y, mask_z, mask_x_loop,  &
              mask_y_loop, mask_z_loop, netcdf_data_format, netcdf_deflate,     &
@@ -719,5 +723,6 @@
              dt_domask, dt_dopr, dt_dopr_listing, dt_dots, dt_do2d_xy,         &
              dt_do2d_xz, dt_do2d_yz, dt_do3d, dt_max, dt_restart,              &
-             dt_run_control,end_time, force_print_header, mask_scale_x,        &
+             dt_run_control,end_time, force_print_header, mask_k_over_surface, &
+             mask_scale_x,        &
              mask_scale_y, mask_scale_z, mask_x, mask_y, mask_z, mask_x_loop,  &
              mask_y_loop, mask_z_loop, netcdf_data_format, netcdf_deflate,     &
Index: /palm/trunk/SOURCE/radiation_model_mod.f90
===================================================================
--- /palm/trunk/SOURCE/radiation_model_mod.f90	(revision 3434)
+++ /palm/trunk/SOURCE/radiation_model_mod.f90	(revision 3435)
@@ -28,4 +28,10 @@
 ! -----------------
 ! $Id$
+! - workaround: return unit=illegal in check_data_output for certain variables
+!   when check called from init_masks
+! - Use pointer in masked output to reduce code redundancies
+! - Add terrain-following masked output
+! 
+! 3424 2018-10-25 07:29:10Z gronemeier
 ! bugfix: add rad_lw_in, rad_lw_out, rad_sw_out to radiation_check_data_output
 ! 
@@ -1250,4 +1256,9 @@
                  'rrtm_asdir*', 'rad_lw_in*', 'rad_lw_out*', 'rad_sw_in*',     &
                  'rad_sw_out*')
+             IF ( i == 0 .AND. ilen == 0 .AND. k == 0)  THEN
+                ! Workaround for masked output (calls with i=ilen=k=0)
+                unit = 'illegal'
+                RETURN
+             ENDIF
              IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
                 message_string = 'illegal value for data_output: "' //         &
@@ -1281,4 +1292,11 @@
 
           CASE ( 'rad_mrt', 'rad_mrt_sw', 'rad_mrt_lw'  )
+
+             IF ( i == 0 .AND. ilen == 0 .AND. k == 0)  THEN
+                ! Workaround for masked output (calls with i=ilen=k=0)
+                unit = 'illegal'
+                RETURN
+             ENDIF
+
              IF ( .NOT.  radiation ) THEN
                 message_string = 'output of "' // TRIM( var ) // '" require'&
@@ -8890,5 +8908,4 @@
     found  = .TRUE.
 
-
 !
 !-- Check for the grid
@@ -9634,10 +9651,15 @@
     CHARACTER (LEN=*) ::  variable   !< 
 
-    INTEGER(iwp) ::  av   !< 
-    INTEGER(iwp) ::  i    !< 
-    INTEGER(iwp) ::  j    !< 
-    INTEGER(iwp) ::  k    !< 
-
-    LOGICAL ::  found     !< 
+    CHARACTER(LEN=5) ::  grid        !< flag to distinquish between staggered grids
+
+    INTEGER(iwp) ::  av              !< 
+    INTEGER(iwp) ::  i               !< 
+    INTEGER(iwp) ::  j               !< 
+    INTEGER(iwp) ::  k               !< 
+    INTEGER(iwp) ::  topo_top_ind    !< k index of highest horizontal surface
+
+    LOGICAL ::  found                !< true if output array was found
+    LOGICAL ::  resorted             !< true if array is resorted
+
 
     REAL(wp),                                                                  &
@@ -9645,6 +9667,10 @@
           local_pf   !< 
 
-
-    found = .TRUE.
+    REAL(wp), DIMENSION(:,:,:), POINTER ::  to_be_resorted  !< points to array which needs to be resorted for output
+
+
+    found    = .TRUE.
+    grid     = 's'
+    resorted = .FALSE.
 
     SELECT CASE ( TRIM( variable ) )
@@ -9653,174 +9679,103 @@
        CASE ( 'rad_lw_in' )
           IF ( av == 0 )  THEN
-             DO  i = 1, mask_size_l(mid,1)
-                DO  j = 1, mask_size_l(mid,2)
-                   DO  k = 1, mask_size_l(mid,3)
-                       local_pf(i,j,k) = rad_lw_in(mask_k(mid,k),              &
-                                            mask_j(mid,j),mask_i(mid,i))
-                    ENDDO
-                 ENDDO
-              ENDDO
+             to_be_resorted => rad_lw_in
           ELSE
-             DO  i = 1, mask_size_l(mid,1)
-                DO  j = 1, mask_size_l(mid,2)
-                   DO  k = 1, mask_size_l(mid,3)
-                       local_pf(i,j,k) = rad_lw_in_av(mask_k(mid,k),           &
-                                               mask_j(mid,j),mask_i(mid,i))
-                   ENDDO
+             to_be_resorted => rad_lw_in_av
+          ENDIF
+
+       CASE ( 'rad_lw_out' )
+          IF ( av == 0 )  THEN
+             to_be_resorted => rad_lw_out
+          ELSE
+             to_be_resorted => rad_lw_out_av
+          ENDIF
+
+       CASE ( 'rad_lw_cs_hr' )
+          IF ( av == 0 )  THEN
+             to_be_resorted => rad_lw_cs_hr
+          ELSE
+             to_be_resorted => rad_lw_cs_hr_av
+          ENDIF
+
+       CASE ( 'rad_lw_hr' )
+          IF ( av == 0 )  THEN
+             to_be_resorted => rad_lw_hr
+          ELSE
+             to_be_resorted => rad_lw_hr_av
+          ENDIF
+
+       CASE ( 'rad_sw_in' )
+          IF ( av == 0 )  THEN
+             to_be_resorted => rad_sw_in
+          ELSE
+             to_be_resorted => rad_sw_in_av
+          ENDIF
+
+       CASE ( 'rad_sw_out' )
+          IF ( av == 0 )  THEN
+             to_be_resorted => rad_sw_out
+          ELSE
+             to_be_resorted => rad_sw_out_av
+          ENDIF
+
+       CASE ( 'rad_sw_cs_hr' )
+          IF ( av == 0 )  THEN
+             to_be_resorted => rad_sw_cs_hr
+          ELSE
+             to_be_resorted => rad_sw_cs_hr_av
+          ENDIF
+
+       CASE ( 'rad_sw_hr' )
+          IF ( av == 0 )  THEN
+             to_be_resorted => rad_sw_hr
+          ELSE
+             to_be_resorted => rad_sw_hr_av
+          ENDIF
+
+       CASE DEFAULT
+          found = .FALSE.
+
+    END SELECT
+
+!
+!-- Resort the array to be output, if not done above
+    IF ( .NOT. resorted )  THEN
+       IF ( .NOT. mask_surface(mid) )  THEN
+!
+!--       Default masked output
+          DO  i = 1, mask_size_l(mid,1)
+             DO  j = 1, mask_size_l(mid,2)
+                DO  k = 1, mask_size_l(mid,3)
+                   local_pf(i,j,k) =  to_be_resorted(mask_k(mid,k), &
+                                      mask_j(mid,j),mask_i(mid,i))
                 ENDDO
              ENDDO
-          ENDIF
-
-       CASE ( 'rad_lw_out' )
-          IF ( av == 0 )  THEN
-             DO  i = 1, mask_size_l(mid,1)
-                DO  j = 1, mask_size_l(mid,2)
-                   DO  k = 1, mask_size_l(mid,3)
-                       local_pf(i,j,k) = rad_lw_out(mask_k(mid,k),             &
-                                            mask_j(mid,j),mask_i(mid,i))
-                    ENDDO
-                 ENDDO
-              ENDDO
-          ELSE
-             DO  i = 1, mask_size_l(mid,1)
-                DO  j = 1, mask_size_l(mid,2)
-                   DO  k = 1, mask_size_l(mid,3)
-                       local_pf(i,j,k) = rad_lw_out_av(mask_k(mid,k),          &
-                                               mask_j(mid,j),mask_i(mid,i))
-                   ENDDO
+          ENDDO
+
+       ELSE
+!
+!--       Terrain-following masked output
+          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
+                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)                     )
                 ENDDO
              ENDDO
-          ENDIF
-
-       CASE ( 'rad_lw_cs_hr' )
-          IF ( av == 0 )  THEN
-             DO  i = 1, mask_size_l(mid,1)
-                DO  j = 1, mask_size_l(mid,2)
-                   DO  k = 1, mask_size_l(mid,3)
-                       local_pf(i,j,k) = rad_lw_cs_hr(mask_k(mid,k),           &
-                                            mask_j(mid,j),mask_i(mid,i))
-                    ENDDO
-                 ENDDO
-              ENDDO
-          ELSE
-             DO  i = 1, mask_size_l(mid,1)
-                DO  j = 1, mask_size_l(mid,2)
-                   DO  k = 1, mask_size_l(mid,3)
-                       local_pf(i,j,k) = rad_lw_cs_hr_av(mask_k(mid,k),        &
-                                               mask_j(mid,j),mask_i(mid,i))
-                   ENDDO
-                ENDDO
-             ENDDO
-          ENDIF
-
-       CASE ( 'rad_lw_hr' )
-          IF ( av == 0 )  THEN
-             DO  i = 1, mask_size_l(mid,1)
-                DO  j = 1, mask_size_l(mid,2)
-                   DO  k = 1, mask_size_l(mid,3)
-                       local_pf(i,j,k) = rad_lw_hr(mask_k(mid,k),              &
-                                            mask_j(mid,j),mask_i(mid,i))
-                    ENDDO
-                 ENDDO
-              ENDDO
-          ELSE
-             DO  i = 1, mask_size_l(mid,1)
-                DO  j = 1, mask_size_l(mid,2)
-                   DO  k = 1, mask_size_l(mid,3)
-                       local_pf(i,j,k) = rad_lw_hr_av(mask_k(mid,k),           &
-                                               mask_j(mid,j),mask_i(mid,i))
-                   ENDDO
-                ENDDO
-             ENDDO
-          ENDIF
-
-       CASE ( 'rad_sw_in' )
-          IF ( av == 0 )  THEN
-             DO  i = 1, mask_size_l(mid,1)
-                DO  j = 1, mask_size_l(mid,2)
-                   DO  k = 1, mask_size_l(mid,3)
-                       local_pf(i,j,k) = rad_sw_in(mask_k(mid,k),              &
-                                            mask_j(mid,j),mask_i(mid,i))
-                    ENDDO
-                 ENDDO
-              ENDDO
-          ELSE
-             DO  i = 1, mask_size_l(mid,1)
-                DO  j = 1, mask_size_l(mid,2)
-                   DO  k = 1, mask_size_l(mid,3)
-                       local_pf(i,j,k) = rad_sw_in_av(mask_k(mid,k),           &
-                                               mask_j(mid,j),mask_i(mid,i))
-                   ENDDO
-                ENDDO
-             ENDDO
-          ENDIF
-
-       CASE ( 'rad_sw_out' )
-          IF ( av == 0 )  THEN
-             DO  i = 1, mask_size_l(mid,1)
-                DO  j = 1, mask_size_l(mid,2)
-                   DO  k = 1, mask_size_l(mid,3)
-                       local_pf(i,j,k) = rad_sw_out(mask_k(mid,k),             &
-                                            mask_j(mid,j),mask_i(mid,i))
-                    ENDDO
-                 ENDDO
-              ENDDO
-          ELSE
-             DO  i = 1, mask_size_l(mid,1)
-                DO  j = 1, mask_size_l(mid,2)
-                   DO  k = 1, mask_size_l(mid,3)
-                       local_pf(i,j,k) = rad_sw_out_av(mask_k(mid,k),          &
-                                               mask_j(mid,j),mask_i(mid,i))
-                   ENDDO
-                ENDDO
-             ENDDO
-          ENDIF
-
-       CASE ( 'rad_sw_cs_hr' )
-          IF ( av == 0 )  THEN
-             DO  i = 1, mask_size_l(mid,1)
-                DO  j = 1, mask_size_l(mid,2)
-                   DO  k = 1, mask_size_l(mid,3)
-                       local_pf(i,j,k) = rad_sw_cs_hr(mask_k(mid,k),           &
-                                            mask_j(mid,j),mask_i(mid,i))
-                    ENDDO
-                 ENDDO
-              ENDDO
-          ELSE
-             DO  i = 1, mask_size_l(mid,1)
-                DO  j = 1, mask_size_l(mid,2)
-                   DO  k = 1, mask_size_l(mid,3)
-                       local_pf(i,j,k) = rad_sw_cs_hr_av(mask_k(mid,k),        &
-                                               mask_j(mid,j),mask_i(mid,i))
-                   ENDDO
-                ENDDO
-             ENDDO
-          ENDIF
-
-       CASE ( 'rad_sw_hr' )
-          IF ( av == 0 )  THEN
-             DO  i = 1, mask_size_l(mid,1)
-                DO  j = 1, mask_size_l(mid,2)
-                   DO  k = 1, mask_size_l(mid,3)
-                       local_pf(i,j,k) = rad_sw_hr(mask_k(mid,k),              &
-                                            mask_j(mid,j),mask_i(mid,i))
-                    ENDDO
-                 ENDDO
-              ENDDO
-          ELSE
-             DO  i = 1, mask_size_l(mid,1)
-                DO  j = 1, mask_size_l(mid,2)
-                   DO  k = 1, mask_size_l(mid,3)
-                       local_pf(i,j,k) = rad_sw_hr_av(mask_k(mid,k),           &
-                                               mask_j(mid,j),mask_i(mid,i))
-                   ENDDO
-                ENDDO
-             ENDDO
-          ENDIF
-
-       CASE DEFAULT
-          found = .FALSE.
-
-    END SELECT
+          ENDDO
+
+       ENDIF
+    ENDIF
+
 
 
Index: /palm/trunk/SOURCE/urban_surface_mod.f90
===================================================================
--- /palm/trunk/SOURCE/urban_surface_mod.f90	(revision 3434)
+++ /palm/trunk/SOURCE/urban_surface_mod.f90	(revision 3435)
@@ -28,4 +28,7 @@
 ! -----------------
 ! $Id$
+! Bugfix: allocate gamma_w_green_sat until nzt_wall+1
+! 
+! 3418 2018-10-24 16:07:39Z kanani
 ! (rvtils, srissman)
 ! -Updated building databse, two green roof types (ind_green_type_roof)
@@ -1376,5 +1379,5 @@
         ALLOCATE ( surf_usm_h%alpha_vg_green(1:surf_usm_h%ns)     )
         ALLOCATE ( surf_usm_h%l_vg_green(1:surf_usm_h%ns)     )
-        ALLOCATE ( surf_usm_h%gamma_w_green_sat(nzb_wall:nzt_wall,1:surf_usm_h%ns)     )
+        ALLOCATE ( surf_usm_h%gamma_w_green_sat(nzb_wall:nzt_wall+1,1:surf_usm_h%ns)     )
         ALLOCATE ( surf_usm_h%lambda_w_green(nzb_wall:nzt_wall,1:surf_usm_h%ns)     )
         ALLOCATE ( surf_usm_h%gamma_w_green(nzb_wall:nzt_wall,1:surf_usm_h%ns)     )
Index: /palm/trunk/SOURCE/user_data_output_mask.f90
===================================================================
--- /palm/trunk/SOURCE/user_data_output_mask.f90	(revision 3434)
+++ /palm/trunk/SOURCE/user_data_output_mask.f90	(revision 3435)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Add terrain-following output
+! 
+! 2718 2018-01-02 08:49:38Z maronga
 ! Corrected "Former revisions" section
 ! 
@@ -62,16 +65,21 @@
     USE kinds
     
+    USE surface_mod,                                                           &
+        ONLY:  get_topography_top_index_ji
+    
     USE user
 
     IMPLICIT NONE
 
-    CHARACTER (LEN=*) ::  variable   !< 
+    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) ::  av             !< 
+    INTEGER(iwp) ::  i              !< 
+    INTEGER(iwp) ::  j              !< 
+    INTEGER(iwp) ::  k              !< 
+    INTEGER(iwp) ::  topo_top_ind   !< k index of highest horizontal surface
 
-    LOGICAL ::  found     !< 
+    LOGICAL ::  found               !< 
 
     REAL(wp),                                                                  &
@@ -81,4 +89,5 @@
 
     found = .TRUE.
+    grid  = 's'
 
     SELECT CASE ( TRIM( variable ) )
@@ -90,21 +99,74 @@
 !       CASE ( 'u2' )
 !          IF ( av == 0 )  THEN
-!            DO  i = 1, mask_size_l(mid,1)
-!               DO  j = 1, mask_size_l(mid,2)
-!                  DO  k = 1, mask_size_l(mid,3)
-!                      local_pf(i,j,k) = u2(mask_k(mid,k),                       &
-!                                           mask_j(mid,j),mask_i(mid,i))
+!             IF ( .NOT. mask_surface(mid) )  THEN
+!!
+!!--             Default masked output
+!                DO  i = 1, mask_size_l(mid,1)
+!                   DO  j = 1, mask_size_l(mid,2)
+!                      DO  k = 1, mask_size_l(mid,3)
+!                         local_pf(i,j,k) = u2(mask_k(mid,k),                  &
+!                                              mask_j(mid,j),                  &
+!                                              mask_i(mid,i))
+!                      ENDDO
 !                   ENDDO
 !                ENDDO
-!             ENDDO
-!          ELSE
-!            DO  i = 1, mask_size_l(mid,1)
-!               DO  j = 1, mask_size_l(mid,2)
-!                  DO  k = 1, mask_size_l(mid,3)
-!                      local_pf(i,j,k) = u2_av(mask_k(mid,k),                    &
-!                                              mask_j(mid,j),mask_i(mid,i))
+!             ELSE
+!!
+!!--             Terrain-following masked output
+!                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
+!                      DO  k = 1, mask_size_l(mid,3)
+!                         local_pf(i,j,k) = u2(MIN( topo_top_ind+mask_k(mid,k),&
+!                                                   nzt+1 ),                   &
+!                                              mask_j(mid,j),                  &
+!                                              mask_i(mid,i)                   )
+!                      ENDDO
 !                   ENDDO
 !                ENDDO
-!             ENDDO
+!             ENDIF
+!          ELSE
+!             IF ( .NOT. mask_surface(mid) )  THEN
+!!
+!!--             Default masked output
+!                DO  i = 1, mask_size_l(mid,1)
+!                   DO  j = 1, mask_size_l(mid,2)
+!                      DO  k = 1, mask_size_l(mid,3)
+!                          local_pf(i,j,k) = u2_av(mask_k(mid,k),              &
+!                                                  mask_j(mid,j),              &
+!                                                  mask_i(mid,i) )
+!                       ENDDO
+!                    ENDDO
+!                 ENDDO
+!             ELSE
+!!
+!!--             Terrain-following masked output
+!                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
+!                      DO  k = 1, mask_size_l(mid,3)
+!                         local_pf(i,j,k) = u2_av(                             &
+!                                              MIN( topo_top_ind+mask_k(mid,k),&
+!                                                   nzt+1 ),                   &
+!                                              mask_j(mid,j),                  &
+!                                              mask_i(mid,i)                   )
+!                      ENDDO
+!                   ENDDO
+!                ENDDO
+!             ENDIF
 !          ENDIF
 
