Index: palm/trunk/SOURCE/data_output_3d.f90
===================================================================
--- palm/trunk/SOURCE/data_output_3d.f90	(revision 2512)
+++ palm/trunk/SOURCE/data_output_3d.f90	(revision 2696)
@@ -1,5 +1,5 @@
 !> @file data_output_3d.f90
 !------------------------------------------------------------------------------!
-! This file is part of PALM.
+! This file is part of the PALM model system.
 !
 ! PALM is free software: you can redistribute it and/or modify it under the
@@ -25,4 +25,10 @@
 ! -----------------
 ! $Id$
+! Implementation of turbulence_closure_mod (TG)
+! Implementation of chemistry module (FK)
+! Set fill values at topography grid points or e.g. non-natural-type surface
+! in case of LSM output (MS)
+! 
+! 2512 2017-10-04 08:26:59Z raasch
 ! upper bounds of 3d output changed from nx+1,ny+1 to nx,ny
 ! no output of ghost layer data
@@ -169,10 +175,16 @@
     USE averaging
         
+#if defined( __chem )
+    USE chemistry_model_mod,                                                   &
+        ONLY:  chem_data_output_3d
+#endif        
+
     USE cloud_parameters,                                                      &
         ONLY:  l_d_cp, pt_d_t
         
     USE control_parameters,                                                    &
-        ONLY:  cloud_physics, do3d, do3d_no, do3d_time_count, io_blocks,       &
-               io_group, land_surface, message_string, ntdim_3d, nz_do3d,      &
+        ONLY:  air_chemistry, cloud_physics, do3d, do3d_no, do3d_time_count,   &
+               io_blocks, io_group, land_surface, message_string,              &
+               ntdim_3d, nz_do3d,                                              &
                psolver, simulated_time, time_since_reference_point,            &
                urban_surface, varnamelength
@@ -182,5 +194,5 @@
         
     USE indices,                                                               &
-        ONLY:  nbgp, nx, nxl, nxr, ny, nyn, nys, nzb
+        ONLY:  nbgp, nx, nxl, nxr, ny, nyn, nys, nzb, wall_flags_0
         
     USE kinds
@@ -194,5 +206,5 @@
 
     USE netcdf_interface,                                                      &
-        ONLY:  id_set_3d, id_var_do3d, id_var_time_3d, nc_stat,                &
+        ONLY:  fill_value, id_set_3d, id_var_do3d, id_var_time_3d, nc_stat,    &
                netcdf_data_format, netcdf_handle_error
         
@@ -207,8 +219,11 @@
         
     USE radiation_model_mod,                                                   &
-        ONLY:  radiation, radiation_data_output_3d
+        ONLY:  nzub, nzut, radiation, radiation_data_output_3d
+
+    USE turbulence_closure_mod,                                                &
+        ONLY:  tcm_data_output_3d
 
     USE urban_surface_mod,                                                     &
-        ONLY:  nzub, nzut, usm_data_output_3d
+        ONLY:  usm_data_output_3d
 
 
@@ -216,4 +231,5 @@
 
     INTEGER(iwp) ::  av        !< 
+    INTEGER(iwp) ::  flag_nr   !< number of masking flag
     INTEGER(iwp) ::  i         !< 
     INTEGER(iwp) ::  if        !< 
@@ -315,4 +331,10 @@
 !--    Allocate a temporary array with the desired output dimensions.
        ALLOCATE( local_pf(nxl:nxr,nys:nyn,nzb_do:nzt_do) )
+!
+!--    Before each output, set array local_pf to fill value
+       local_pf = fill_value
+!
+!--    Set masking flag for topography for not resorted arrays
+       flag_nr = 0
 
        SELECT CASE ( trimvar )
@@ -569,4 +591,5 @@
 
           CASE ( 'u' )
+             flag_nr = 1
              IF ( av == 0 )  THEN
                 to_be_resorted => u
@@ -576,4 +599,5 @@
 
           CASE ( 'v' )
+             flag_nr = 2
              IF ( av == 0 )  THEN
                 to_be_resorted => v
@@ -590,4 +614,5 @@
 
           CASE ( 'w' )
+             flag_nr = 3
              IF ( av == 0 )  THEN
                 to_be_resorted => w
@@ -613,4 +638,5 @@
                 DEALLOCATE ( local_pf )
                 ALLOCATE( local_pf(nxl:nxr,nys:nyn,nzb_do:nzt_do) )
+                local_pf = fill_value
 
                 CALL lsm_data_output_3d( av, do3d(av,if), found, local_pf )
@@ -629,4 +655,9 @@
              ENDIF
 
+             IF ( .NOT. found )  THEN
+                CALL tcm_data_output_3d( av, do3d(av,if), found, local_pf )
+                resorted = .TRUE.
+             ENDIF
+
 !
 !--          Radiation quantity
@@ -638,8 +669,18 @@
 
 !
+!--          Chemistry model output
+#if defined( __chem )
+             IF ( .NOT. found  .AND.  air_chemistry )  THEN
+                CALL chem_data_output_3d( av, do3d(av,if), found,              &
+                                          local_pf )
+                resorted = .TRUE.
+             ENDIF
+#endif
+
+!
 !--          Plant canopy model output
              IF ( .NOT. found  .AND.  plant_canopy )  THEN
-                CALL pcm_data_output_3d( av, do3d(av,if), found,         &
-                                               local_pf )
+                CALL pcm_data_output_3d( av, do3d(av,if), found, local_pf,     &
+                                         fill_value )
                 resorted = .TRUE.
              ENDIF
@@ -667,5 +708,8 @@
              DO  j = nys, nyn
                 DO  k = nzb_do, nzt_do
-                   local_pf(i,j,k) = to_be_resorted(k,j,i)
+                   local_pf(i,j,k) = MERGE(                                    &
+                                      to_be_resorted(k,j,i),                   &
+                                      REAL( fill_value, KIND = wp ),           &
+                                      BTEST( wall_flags_0(k,j,i), flag_nr ) )
                 ENDDO
              ENDDO
