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