Index: /palm/trunk/SCRIPTS/create_basic_static_driver.py
===================================================================
--- /palm/trunk/SCRIPTS/create_basic_static_driver.py	(revision 4795)
+++ /palm/trunk/SCRIPTS/create_basic_static_driver.py	(revision 4796)
@@ -1,316 +1,458 @@
-#------------------------------------------------------------------------------#
+#!/usr/bin/env python3
+# -------------------------------------------------------------------- #
 # This file is part of the PALM model system.
 #
-# PALM is free software: you can redistribute it and/or modify it under the
-# terms of the GNU General Public License as published by the Free Software
-# Foundation, either version 3 of the License, or (at your option) any later
-# version.
+# PALM is free software: you can redistribute it and/or modify it under
+# the terms of the GNU General Public License as published by the Free
+# Software Foundation, either version 3 of the License, or (at your
+# option) any later version.
 #
-# PALM is distributed in the hope that it will be useful, but WITHOUT ANY
-# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
-# A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
+# PALM is distributed in the hope that it will be useful, but WITHOUT
+# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+# FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+# for more details.
 #
-# You should have received a copy of the GNU General Public License along with
-# PALM. If not, see <http://www.gnu.org/licenses/>.
+# You should have received a copy of the GNU General Public License
+# along with PALM. If not, see <http://www.gnu.org/licenses/>.
 #
 # Copyright 1997-2020 Leibniz Universitaet Hannover
-#------------------------------------------------------------------------------#
-
+# -------------------------------------------------------------------- #
+
+import datetime
 import math
-from   netCDF4 import Dataset
+import pytz
+
+from netCDF4 import Dataset
 import numpy as np
-import datetime
 
 
 class StaticDriver:
-    """ This is an example script to generate static drivers for PALM. 
-    
+    """This is an example script to generate static drivers for PALM.
+
     You can use it as a starting point for creating your setup specific
-    driver.  
+    driver.
     """
-    
-    
+
     def __init__(self):
-        """ Open the static driver as NetCDF4 file. Here, you have to give the 
-        full path to the static driver that shall be created. Existing file 
-        with same name is deleted. 
+        """Open the static driver as netCDF4 file. Here, you have to
+        give the full path to the static driver that shall be created.
+        Existing file with same name is deleted.
         """
         print('Opening file...')
-        self.nc_file = Dataset('path/to/file.nc', 'w', format='NETCDF4')
-
+        self.nc_file = Dataset('example_static_file.nc', 'w', format='NETCDF4')
 
     def write_global_attributes(self):
-        """ Write global attributes to static driver. """
+        """Write global attributes to static driver."""
         print("Writing global attributes...")
-        
-        # mandatory global attributes
-        self.nc_file.origin_lon = 55.0 # used to initialize coriolis parameter
-        self.nc_file.origin_lat = 0.0 # (overwrite initialization_parameters)
-        self.nc_file.origin_time = '2000-06-21 12:00:00 +00'
-        self.nc_file.origin_x = 308124
-        self.nc_file.origin_y = 6098908
+
+        # Optional global attributes
+        # --------------------------
+        self.nc_file.title = 'Example PALM static driver'
+        self.nc_file.author = 'PALM user'
+        self.nc_file.institution = 'Institut of Meteorology and Climatology,' \
+            'Leibniz University Hannover'
+        self.nc_file.comment = 'Generic crossing example'
+        self.nc_file.creation_date = \
+            pytz.utc.localize(datetime.datetime.utcnow()).strftime('%y-%m-%d %H:%M:%S %z')[:-2]
+        self.nc_file.history = ''
+        self.nc_file.keywords = 'example, PALM-4U'
+        self.nc_file.license = ''
+        self.nc_file.palm_version = ''
+        self.nc_file.references = ''
+        self.nc_file.source = 'PALM trunk'
+        self.nc_file.version = '1'
+
+        # Mandatory global attributes
+        # ---------------------------
+        self.nc_file.Conventions = 'CF-1.7'
+        self.nc_file.origin_lat = 52.50965  # (overwrite initialization_parameters)
+        self.nc_file.origin_lon = 13.3139  # Used to initialize Coriolis parameter
+        self.nc_file.origin_time = '2019-03-06 10:00:00 +00'
+        self.nc_file.origin_x = 3455249.0
+        self.nc_file.origin_y = 5424815.0
         self.nc_file.origin_z = 0.0
         self.nc_file.rotation_angle = 0.0
-        
-        # optional global attributes
-        self.nc_file.author = 'Your Name'
-        self.nc_file.comment = 'Miscellaneous information about the data ' \
-                               'or methods to produce it.'
-        self.nc_file.creation_date = str(datetime.datetime.now())
-        self.nc_file.institution = 'Institut of Meteorology and Climatology,' \
-                                   'Leibniz University Hannover'
-        self.nc_file.history = ''
-        self.nc_file.palm_revision = ''
-        self.nc_file.title = 'Static driver for some arbitrary PALM setup'
-
 
     def define_dimensions(self):
-        """ Set dimensions on which variables are defined. """
+        """Set dimensions on which variables are defined."""
         print("Writing dimensions...")
-        
-        # specify general grid parameters
-        # these values must equal those set in the initialization_parameters
-        self.nx = 39
-        self.ny = 39
-        self.nz = 40
-        dx = 50
-        dy = 50
-        dz = 50
-        
-        # create soil grid (only relevant if land surface module is used)
+
+        # Specify general grid parameters
+        # These values must equal to those set in the initialization_parameters
+        self.nx = 19
+        self.ny = 19
+        self.nz = 20
+        dx = 2
+        dy = 2
+        dz = 2
+
+        # Create soil grid (only relevant if land surface module is used)
         dz_soil = np.array((0.01, 0.02, 0.04, 0.06, 0.14, 0.26, 0.54, 1.86))
         zsoil_fullLayers = np.zeros_like(dz_soil)
-        zsoil_fullLayers = np.around([np.sum(dz_soil[:zs]) for zs in np.arange(1,len(dz_soil)+1)],2)
+        zsoil_fullLayers = np.around(
+            [np.sum(dz_soil[:zs]) for zs in np.arange(1, len(dz_soil)+1)], 2)
         zsoil_array = zsoil_fullLayers - dz_soil/2.
-        
-        # create plant canopy grid (only relevant of plant canopy module is used)
-        zlad_array = np.arange(0,10,1)
-        
-        # mandatory dimensions
-        self.nc_file.createDimension('x' ,self.nx+1)
-        self.x = self.nc_file.createVariable('x', 'f8', ('x',))
+
+        # Coordinates
+        # -----------
+        self.nc_file.createDimension('x', self.nx+1)
+        self.x = self.nc_file.createVariable('x', 'f4', ('x',))
+        self.x.long_name = 'distance to origin in x-direction'
         self.x.units = 'm'
-        self.x.standard_name = 'x coordinate of cell centers'
-        self.x[:] = np.arange(0,(self.nx+1)*dx,dx)
-        
-        self.nc_file.createDimension('xu', self.nx)
-        self.xu = self.nc_file.createVariable('xu', 'f8', ('xu',))
-        self.xu.units = 'm'
-        self.xu.standard_name = 'x coordinate of cell edges'
-        self.xu[:] = np.arange(dx/2.,self.nx*dx,dx)
-        
+        self.x.axis = 'X'
+        self.x[:] = np.arange(0, (self.nx+1)*dx, dx) + 0.5 * dx
+
         self.nc_file.createDimension('y', self.ny+1)
-        self.y = self.nc_file.createVariable('y', 'f8', ('y',))
+        self.y = self.nc_file.createVariable('y', 'f4', ('y',))
+        self.y.long_name = 'distance to origin in y-direction'
         self.y.units = 'm'
-        self.y.standard_name = 'y coordinate of cell centers'
-        self.y[:] = np.arange(0,(self.ny+1)*dy,dy)
-        
-        self.nc_file.createDimension('yv', self.ny)
-        self.yv = self.nc_file.createVariable('yv', 'f8', ('yv',))
-        self.yv.units = 'm'
-        self.yv.standard_name = 'y coordinate of cell edges'
-        self.yv[:] = np.arange(dy/2.,self.ny*dy,dy)
-        
-        # if your simulation uses a stretched vertical grid, you need to 
-        # modify the z and zw coordinates, 
-        # e.g. z_array = (...) and zw_array = (...)
-        z_array = np.append(0, np.arange(dz/2,(self.nz+1)*dz,dz))
-        self.nc_file.createDimension('z',self.nz+2)
-        self.z = self.nc_file.createVariable('z', 'f8', ('z',))
+        self.y.axis = 'Y'
+        self.y[:] = np.arange(0, (self.ny+1)*dy, dy) + 0.5 * dy
+
+        # NOTE if your simulation uses a stretched vertical grid, you
+        # need to modify the z coordinates, e.g. z_array = (...)
+        z_array = np.append(0, np.arange(dz/2, (self.nz)*dz, dz))
+        self.nc_file.createDimension('z', self.nz+1)
+        self.z = self.nc_file.createVariable('z', 'f4', ('z',))
+        self.z.long_name = 'height above origin'
         self.z.units = 'm'
-        self.z.standard_name = 'z coordinate of cell centers'
+        self.z.axis = 'Z'
+        self.z.positive = 'up'
         self.z[:] = z_array
 
-        zw_array = np.arange(0,(self.nz+2)*dz,dz)
-        self.nc_file.createDimension('zw', self.nz+2)
-        self.zw = self.nc_file.createVariable('zw', 'f8', ('zw',))
-        self.zw.units = 'm'
-        self.zw.standard_name = 'z coordinate of cell edges'
-        self.zw[:] = zw_array
-        
-        # optional dimensions, uncomment if needed
-        #self.nc_file.createDimension('zsoil', len(zsoil_array))
-        #self.zsoil = self.nc_file.createVariable('zsoil', 'f8', ('zsoil',))
-        #self.zsoil.positive = 'down'
-        #self.zsoil.units = 'm'
-        #self.zsoil.standard_name = 'depth below land'
-        #self.zsoil[:] = zsoil_array
-        
-        #self.nc_file.createDimension('zlad', len(zlad_array))
-        #self.zlad = self.nc_file.createVariable('zlad', 'f8', ('zlad',))
-        #self.zlad.units = 'm'
-        #self.zlad.standard_name = 'z coordinate of resolved plant canopy'
-        #self.zlad[:] = zlad_array
-        
-        #self.nc_file.createDimension('nalbedo_pars', 3)
-        #self.nalbedo_pars = self.nc_file.createVariable(
-            #'nalbedo_pars', 'i1', ('nalbedo_pars',))
-        #self.nalbedo_pars[:] = np.arange(0,3)
-        
-        #self.nc_file.createDimension('nbuilding_pars', 46)
-        #self.nbuilding_pars = self.nc_file.createVariable(
-            #'nbuilding_pars', 'i4', ('nbuilding_pars',))
-        #self.nbuilding_pars[:] = np.arange(0,46)
-        
-        #self.nc_file.createDimension('nsoil_pars', 8)
-        #self.nsoil_pars = self.nc_file.createVariable(
-            #'nsoil_pars', 'i1', ('nsoil_pars',))
-        #self.nsoil_pars[:] = np.arange(0,8)
-        
-        #self.nc_file.createDimension('nsurface_fraction', 3)
-        #self.nsurface_fraction = self.nc_file.createVariable(
-            #'nsurface_fraction', 'i1', ('nsurface_fraction',))
-        #self.nsurface_fraction[:] = np.arange(0,3)
-        
-        #self.nc_file.createDimension('nvegetation_pars', 12)
-        #self.nvegetation_pars = self.nc_file.createVariable(
-            #'nvegetation_pars', 'i1', ('nvegetation_pars',))
-        #self.nvegetation_pars[:] = np.arange(0,12)
-        
-        #self.nc_file.createDimension('npavement_pars', 4)
-        #self.npavement_pars = self.nc_file.createVariable(
-            #'npavement_pars', 'i1', ('npavement_pars',))
-        #self.npavement_pars[:] = np.arange(0,4)
-
-        #self.nc_file.createDimension('npavement_subsurface_pars', 2)
-        #self.npavement_subsurface_pars = self.nc_file.createVariable(
-            #'npavement_subsurface_pars', 'i1', ('npavement_subsurface_pars',))
-        #self.npavement_subsurface_pars[:] = np.arange(0,2)
-        
-        #self.nc_file.createDimension('nwater_pars', 6)
-        #self.nwater_pars = self.nc_file.createVariable(
-            #'nwater_pars', 'i1', ('nwater_pars',))
-        #self.nwater_pars[:] = np.arange(0,6)
-        
-
-    def add_variables(self):
-        """ Uncomment variables below as you like. 
-        
-        Be aware that some variables depend on others. For a description of 
-        each variable, please have a look at the documentation at
+        zlad_array = self.z[:6]
+        self.nc_file.createDimension('zlad', len(zlad_array))
+        self.zlad = self.nc_file.createVariable('zlad', 'f4', ('zlad',))
+        self.zlad.long_name = 'height above ground'
+        self.zlad.units = 'm'
+        self.zlad.axis = 'Z'
+        self.zlad.positive = 'up'
+        self.zlad[:] = zlad_array
+
+        # self.nc_file.createDimension('zsoil', len(zsoil_array))
+        # self.zsoil = self.nc_file.createVariable('zsoil', 'f4', ('zsoil',))
+        # self.zsoil.long_name = 'depth in the soil'
+        # self.zsoil.units = 'm'
+        # self.zsoil.axis = 'Z'
+        # self.zsoil.positive = 'down'
+        # self.zsoil[:] = zsoil_array
+
+        # Other dimensions
+        # ----------------
+        # self.nc_file.createDimension('nalbedo_pars', 3)
+        # self.nalbedo_pars = self.nc_file.createVariable(
+        #     'nalbedo_pars', 'i4', ('nalbedo_pars',))
+        # self.nalbedo_pars[:] = np.arange(8)
+        #
+        # self.nc_file.createDimension('nbuilding_pars', 46)
+        # self.nbuilding_pars = self.nc_file.createVariable(
+        #     'nbuilding_pars', 'i4', ('nbuilding_pars',))
+        # self.nbuilding_pars[:] = np.arange(46)
+        #
+        # self.nc_file.createDimension('npavement_pars', 4)
+        # self.npavement_pars = self.nc_file.createVariable(
+        #     'npavement_pars', 'i4', ('npavement_pars',))
+        # self.npavement_pars[:] = np.arange(4)
+        #
+        # self.nc_file.createDimension('npavement_subsurface_pars', 2)
+        # self.npavement_subsurface_pars = self.nc_file.createVariable(
+        #     'npavement_subsurface_pars', 'i4', ('npavement_subsurface_pars',))
+        # self.npavement_subsurface_pars[:] = np.arange(2)
+        #
+        # self.nc_file.createDimension('nsoil_pars', 8)
+        # self.nsoil_pars = self.nc_file.createVariable(
+        #     'nsoil_pars', 'i4', ('nsoil_pars',))
+        # self.nsoil_pars[:] = np.arange(8)
+        #
+        self.nc_file.createDimension('nsurface_fraction', 3)
+        self.nsurface_fraction = self.nc_file.createVariable(
+            'nsurface_fraction', 'i4', ('nsurface_fraction',))
+        self.nsurface_fraction[:] = np.arange(3)
+
+        # self.nc_file.createDimension('nvegetation_pars', 12)
+        # self.nvegetation_pars = self.nc_file.createVariable(
+        #     'nvegetation_pars', 'i4', ('nvegetation_pars',))
+        # self.nvegetation_pars[:] = np.arange(12)
+        #
+        # self.nc_file.createDimension('nwater_pars', 6)
+        # self.nwater_pars = self.nc_file.createVariable(
+        #     'nwater_pars', 'i4', ('nwater_pars',))
+        # self.nwater_pars[:] = np.arange(7)
+
+    def define_variables(self):
+        """Define variables for the static driver.
+
+        Be aware that some variables depend on others. For a description
+        of each variable, please have a look at the documentation at
         palm.muk.uni-hannover.de/trac/wiki/doc/app/iofiles/pids/static
-        
+
         An example of how you modify the variables is given below:
-        
-        building_2d_array = np.ones((self.ny+1,self.nx+1)) * -9999.9
+
+        building_2d_array = np.ones((self.ny+1, self.nx+1)) * -9999.0
         south_wall, north_wall, left_wall, right_wall = 20, 25, 20, 25
-        building_2d_array[south_wall:north_wall,left_wall:right_wall] = 50
+        building_2d_array[
+            south_wall:north_wall,
+            left_wall:right_wall
+            ] = 50
         nc_buildings_2d = self.nc_file.createVariable(
-            'buildings_2d', 'f4', ('y','x'),fill_value=-9999.9)
+            'buildings_2d', 'f4', ('y', 'x'), fill_value=-9999.0)
         nc_buildings_2d.lod = 1
-        nc_buildings_2d[:,:] = building_2d_array
+        nc_buildings_2d[:, :] = building_2d_array
         """
         print("Writing variables...")
-        
-        ## topography variables
-        #nc_buildings_2d = self.nc_file.createVariable(
-            #'buildings_2d', 'f4', ('y','x'), fill_value=-9999.9)
-        #nc_buildings_2d.lod = 1
-        #nc_buildings_2d[:,:] = -9999.9
-        
-        #nc_buildings_3d = self.nc_file.createVariable(
-            #'buildings_3d', 'i1', ('z','y','x'),fill_value=-127)
-        #nc_buildings_3d.lod = 2 
-        #nc_buildings_3d[:,:,:] = -127
-        
-        #nc_building_id = self.nc_file.createVariable(
-            #'building_id', 'i1', ('y','x'), fill_value=-127)
-        #nc_building_id[:,:] = -127
-        
-        #nc_zt = self.nc_file.createVariable(
-            #'zt', 'f4', ('y','x'), fill_value=-9999.9)
-        #nc_zt.long_name = 'orography'
-        #nc_zt.units = 'm'
-        #nc_zt[:,:] = -9999.9
-        
-        ## surface variables
-        #nc_albedo_pars = self.nc_file.createVariable(
-            #'albedo_pars', 'f4', ('nalbedo_pars','y','x'), fill_value=-9999.9)
-        #nc_albedo_pars.long_name = 'albedo parameters'
-        #nc_albedo_pars[:,:,:] = -9999.9
-        
-        #nc_albedo_type = self.nc_file.createVariable(
-            #'albedo_type', 'f4', ('y','x'), fill_value=-9999.9)
-        #nc_albedo_type[:,:] = -9999.9
-        
-        #nc_building_pars = self.nc_file.createVariable(
-            #'building_pars', 'f4', ('nbuilding_pars','y','x'), fill_value=-9999.9)
-        #nc_building_pars.long_name = 'building parameters'
-        #nc_building_pars[:,:,:] = -9999.9
-        
-        #nc_building_type = self.nc_file.createVariable(
-            #'building_type', 'i1', ('y','x'), fill_value=-127)
-        #nc_building_type[:,:] = -127
-        
-        #nc_pavement_pars = self.nc_file.createVariable(
-            #'pavement_pars', 'f4', ('npavement_pars','y','x'), fill_value=-9999.9)
-        #nc_pavement_pars[:,:,:] = -9999.9
-        
-        #nc_pavement_subsurface_pars = self.nc_file.createVariable(
-            #'pavement_subsurface_pars', 'i1', 
-            #('npavement_subsurface_pars','zsoil','y','x'), fill_value=-9999.9)
-        #nc_pavement_subsurface_pars[:,:,:,:] = -9999.9
-        
-        #nc_pavement_type = self.nc_file.createVariable(
-            #'pavement_type', 'i1', ('y','x'), fill_value=-127)
-        #nc_pavement_type[:,:] = -127
-        
-        #nc_soil_pars = self.nc_file.createVariable(
-            #'soil_pars', 'f', ('nsoil_pars','zsoil','y','x'), fill_value=-9999.9)
-        #nc_soil_pars.long_name = 'soil parameters'
-        #nc_soil_pars.lod = 2 # if set to 1, adjust dimensions of soil_pars
-        #nc_soil_pars[:,:,:,:] = -9999.9
-        
-        #nc_soil_type = self.nc_file.createVariable(
-            #'soil_type', 'i1', ('zsoil','y','x'), fill_value=-127)
-        #nc_soil_type.lod = 2 # if set to 1, adjust dimensions of soil_type
-        #nc_soil_type[:,:,:] = -127
-        
-        ## required, if several surface types are defined at the same location
-        #nc_surface_fraction = self.nc_file.createVariable(
-            #'surface_fraction', 'f', ('nsurface_fraction','y','x'), fill_value=-9999.9)
-        #nc_surface_fraction[0,:,:] = -9999.9 # vegetation fraction 
-        #nc_surface_fraction[1,:,:] = -9999.9 # pavement fraction
-        #nc_surface_fraction[2,:,:] = -9999.9 # water fraction 
-        
-        #nc_vegetation_pars = self.nc_file.createVariable(
-            #'vegetation_pars', 'f', ('nvegetation_pars','y','x'), fill_value=-9999.9)
-        #nc_vegetation_pars.long_name = 'vegetation parameters'
-        #nc_vegetation_pars[:,:,:] = -9999.9
-        
-        #nc_vegetation_type = self.nc_file.createVariable(
-            #'vegetation_type', 'i1', ('y','x'), fill_value=-127)
-        #nc_vegetation_type[:,:] = -127
-        
-        #nc_street_type = self.nc_file.createVariable(
-            #'street_type', 'i1', ('y','x'), 
-            #fill_value=-127)
-        #nc_street_type[:,:] = -127
-        
-        #nc_water_pars = self.nc_file.createVariable(
-            #'water_pars', 'i1', ('nwater_pars','y','x'), fill_value=-9999.9)
-        #nc_water_pars[:,:,:] = -9999.9
-        
-        #nc_water_type = self.nc_file.createVariable(
-            #'water_type', 'i1',('y','x'), fill_value=-127)
-        #nc_water_type[:,:] = -127
-        
-        ## resolved plant canopy
-        #nc_lad = self.nc_file.createVariable(
-            #'lad', 'f4', ('zlad','y','x'), fill_value=-9999.9)
-        #nc_lad.long_name = 'leaf area density'
-        #nc_lad[:,:,:] = -9999.9
-        
-        #nc_root_area_dens_s = self.nc_file.createVariable(
-            #'root_area_dens_s', 'f4', ('zsoil','y','x'), fill_value=-9999.9)
-        #nc_root_area_dens_s.long_name = 'parameterized root area density'
-        #nc_root_area_dens_s.units = '1/m'
-        #nc_root_area_dens_s[:,:,:] = -9999.9
-        
-        
+
+        # Topography set-up
+        # -----------------
+        nc_building_id = self.nc_file.createVariable(
+            'building_id', 'i4', ('y', 'x'), fill_value=-9999)
+        nc_building_id.long_name = "building id number"
+        nc_building_id.units = "1"
+        nc_building_id[:, :] = nc_building_id._FillValue
+
+        nc_buildings_2d = self.nc_file.createVariable(
+            'buildings_2d', 'f4', ('y', 'x'), fill_value=-9999.0)
+        nc_buildings_2d.long_name = "building height"
+        nc_buildings_2d.units = "m"
+        nc_buildings_2d.lod = np.int32(1)
+        nc_buildings_2d[:, :] = nc_buildings_2d._FillValue
+
+        nc_buildings_3d = self.nc_file.createVariable(
+            'buildings_3d', 'i1', ('z', 'y', 'x'), fill_value=-127)
+        nc_buildings_3d.long_name = "building flag"
+        nc_buildings_3d.units = "1"
+        nc_buildings_3d.lod = np.int32(2)
+        nc_buildings_3d[:, :, :] = nc_buildings_3d._FillValue
+
+        nc_zt = self.nc_file.createVariable(
+            'zt', 'f4', ('y', 'x'), fill_value=-9999.0)
+        nc_zt.long_name = 'terrain height'
+        nc_zt.units = 'm'
+        nc_zt[:, :] = nc_zt._FillValue
+
+        # nc_z0 = self.nc_file.createVariable(
+        #     'z0', 'f4', ('y', 'x'), fill_value=-9999.0)
+        # nc_z0.long_name = 'roughness length for momentum'
+        # nc_z0.units = 'm'
+        # nc_z0[:, :] = nc_z0._FillValue
+
+        # Main surface clasification
+        # --------------------------
+        # nc_albedo_type = self.nc_file.createVariable(
+        #     'albedo_type', 'f4', ('y', 'x'), fill_value=-9999.0)
+        # nc_albedo_type.long_name = "albedo type"
+        # nc_albedo_type.units = "1"
+        # nc_albedo_type[:, :] = nc_albedo_type._FillValue
+        #
+        nc_building_type = self.nc_file.createVariable(
+            'building_type', 'i1', ('y', 'x'), fill_value=-127)
+        nc_building_type.long_name = "building type classification"
+        nc_building_type.units = "1"
+        nc_building_type[:, :] = nc_building_type._FillValue
+
+        nc_pavement_type = self.nc_file.createVariable(
+            'pavement_type', 'i1', ('y', 'x'), fill_value=-127)
+        nc_pavement_type.long_name = "pavement type classification"
+        nc_pavement_type.units = "1"
+        nc_pavement_type[:, :] = nc_pavement_type._FillValue
+
+        nc_soil_type = self.nc_file.createVariable(
+            'soil_type', 'i1', ('y', 'x'), fill_value=-127)
+        nc_soil_type.long_name = "soil type classification"
+        nc_soil_type.units = "1"
+        nc_soil_type.lod = np.int32(1)
+        nc_soil_type[:, :] = nc_soil_type._FillValue
+
+        # NOTE alternatively, define different soil types for different
+        # soil layers
+        # nc_soil_type = self.nc_file.createVariable(
+        #     'soil_type', 'i1', ('zsoil', 'y', 'x'), fill_value=-127)
+        # nc_soil_type.long_name = "soil type"
+        # nc_soil_type.units = "1"
+        # nc_soil_type.lod = np.int32(2)
+        # nc_soil_type[:, :, :] = nc_soil_type._FillValue
+        #
+        # nc_street_crossing = self.nc_file.createVariable(
+        #     'street_crossing', 'i1', ('y', 'x'),
+        #     fill_value=-127)
+        # nc_street_crossing.long_name = "street crossing"
+        # nc_street_crossing.units = "1"
+        # nc_street_crossing.valid_range = np.byte(1)
+        # nc_street_crossing[:, :] = nc_street_crossing._FillValue
+        #
+        nc_street_type = self.nc_file.createVariable(
+            'street_type', 'i1', ('y', 'x'),
+            fill_value=-127)
+        nc_street_type.long_name = "street type classification"
+        nc_street_type.units = "1"
+        nc_street_type[:, :] = nc_street_type._FillValue
+
+        nc_surface_fraction = self.nc_file.createVariable(
+            'surface_fraction', 'f4', ('nsurface_fraction', 'y', 'x'), fill_value=-9999.0)
+        nc_surface_fraction.long_name = "surface fraction"
+        nc_surface_fraction.units = "1"
+        nc_surface_fraction[0, :, :] = nc_surface_fraction._FillValue  # vegetation fraction
+        nc_surface_fraction[1, :, :] = nc_surface_fraction._FillValue  # pavement fraction
+        nc_surface_fraction[2, :, :] = nc_surface_fraction._FillValue  # water fraction
+
+        nc_vegetation_type = self.nc_file.createVariable(
+            'vegetation_type', 'i1', ('y', 'x'), fill_value=-127)
+        nc_vegetation_type.long_name = "vegetation type classification"
+        nc_vegetation_type.units = "1"
+        nc_vegetation_type[:, :] = nc_vegetation_type._FillValue
+
+        nc_water_type = self.nc_file.createVariable(
+            'water_type', 'i1', ('y', 'x'), fill_value=-127)
+        nc_water_type.long_name = "water type classification"
+        nc_water_type.units = "1"
+        nc_water_type[:, :] = nc_water_type._FillValue
+
+        # Pixel-based surface and sub-surface parameters
+        # ----------------------------------------------
+        # nc_albedo_pars = self.nc_file.createVariable(
+        #     'albedo_pars', 'f4', ('nalbedo_pars', 'y', 'x'), fill_value=-9999.0)
+        # nc_albedo_pars.long_name = 'albedo parameters'
+        # nc_albedo_pars[:, :, :] = nc_albedo_pars._FillValue
+        #
+        # nc_building_pars = self.nc_file.createVariable(
+        #     'building_pars', 'f4', ('nbuilding_pars', 'y', 'x'), fill_value=-9999.0)
+        # nc_building_pars.long_name = 'building parameters'
+        # nc_building_pars[:, :, :] = nc_building_pars._FillValue
+        #
+        # nc_pavement_pars = self.nc_file.createVariable(
+        #     'pavement_pars', 'f4', ('npavement_pars', 'y', 'x'), fill_value=-9999.0)
+        # nc_pavement_pars.long_name = "pavement parameters"
+        # nc_pavement_pars[:, :, :] = nc_pavement_pars._FillValue
+        #
+        # nc_pavement_subsurface_pars = self.nc_file.createVariable(
+        #     'pavement_subsurface_pars', 'i1',
+        #     ('npavement_subsurface_pars','zsoil', 'y', 'x'), fill_value=-9999.0)
+        # nc_pavement_subsurface_pars.long_name = "pavement subsurface parameters"
+        # nc_pavement_subsurface_pars[:, :, :, :] = nc_pavement_subsurface_pars._FillValue
+        #
+        # nc_soil_pars = self.nc_file.createVariable(
+        #     'soil_pars', 'f', ('nsoil_pars','zsoil', 'y', 'x'), fill_value=-9999.0)
+        # nc_soil_pars.long_name = "soil parameters"
+        # nc_soil_pars.lod = np.int32(2) # if set to 1, adjust dimensions of soil_pars
+        # nc_soil_pars[:, :, :, :] = nc_soil_pars._FillValue
+        #
+        # nc_vegetation_pars = self.nc_file.createVariable(
+        #     'vegetation_pars', 'f', ('nvegetation_pars', 'y', 'x'), fill_value=-9999.0)
+        # nc_vegetation_pars.long_name = 'vegetation parameters'
+        # nc_vegetation_pars[:, :, :] = nc_vegetation_pars._FillValue
+        #
+        # nc_water_pars = self.nc_file.createVariable(
+        #     'water_pars', 'i1', ('nwater_pars', 'y', 'x'), fill_value=-9999.0)
+        # nc_water_pars.long_name = "water parameters"
+        # nc_water_pars[:, :, :] = nc_water_pars._FillValue
+
+        # Vegetation parameters
+        # ---------------------
+        nc_lad = self.nc_file.createVariable(
+            'lad', 'f4', ('zlad', 'y', 'x'), fill_value=-9999.0)
+        nc_lad.long_name = "leaf area density"
+        nc_lad.units = "m2 m-3"
+        nc_lad[:, :, :] = nc_lad._FillValue
+
+        # nc_bad = self.nc_file.createVariable(
+        #     'bad', 'f4', ('zlad', 'y', 'x'), fill_value=-9999.0)
+        # nc_bad.long_name = "basal area density"
+        # nc_bad.units = "m2 m-3"
+        # nc_bad[:, :, :] = nc_bad._FillValue
+        #
+        # nc_root_area_dens_r = self.nc_file.createVariable(
+        #     'root_area_dens_s', 'f4', ('zsoil', 'y', 'x'), fill_value=-9999.0)
+        # nc_root_area_dens_r.long_name = 'root area density of resolved vegetation'
+        # nc_root_area_dens_r.units = '1'
+        # nc_root_area_dens_r[:, :, :] = nc_root_area_dens_r._FillValue
+        #
+        # nc_root_area_dens_s = self.nc_file.createVariable(
+        #     'root_area_dens_s', 'f4', ('zsoil', 'y', 'x'), fill_value=-9999.0)
+        # nc_root_area_dens_s.long_name = 'root area density of parameterized vegetation'
+        # nc_root_area_dens_s.units = '1'
+        # nc_root_area_dens_s[:, :, :] = nc_root_area_dens_s._FillValue
+        #
+        # nc_tree_id = self.nc_file.createVariable(
+        #     'tree_id', 'i4', ('y', 'x'), fill_value=-9999.0)
+        # nc_tree_id.long_name = "tree id"
+        # nc_tree_id.units = "1"
+        # nc_tree_id[:, :, :] = nc_tree_id._FillValue
+
+        # Set topography
+        # --------------
+        nc_building_id[:5, :5] = 1
+        nc_building_id[-5:, :5] = 2
+        nc_building_id[:5, -5:] = 3
+        nc_building_id[-5:, -5:] = 4
+
+        nc_buildings_2d[:, :] = np.where(
+            nc_building_id[:, :] > nc_building_id._FillValue,
+            nc_building_id[:, :] * 10.0,
+            nc_buildings_2d[:, :])
+
+        for k in range(self.nz+1):
+            nc_buildings_3d[k, :, :] = np.where(nc_buildings_2d[:, :] > self.z[k], 1, 0)
+
+        nc_zt[:, :] = 4.0
+
+        # Set surface types
+        # -----------------
+        nc_building_type[:, :] = np.where(
+            nc_building_id[:, :] > nc_building_id._FillValue,
+            nc_building_id[:, :] * 1,
+            nc_building_type[:, :])
+
+        nc_pavement_type[9:11, :5] = 1
+        nc_pavement_type[:, 7:13] = 2
+
+        nc_vegetation_type[:, 5:7] = 3
+        nc_vegetation_type[:, 13:15] = 3
+        nc_vegetation_type[5:15, 15:] = 3
+        nc_vegetation_type[7:9, 15:17] = 1
+        nc_vegetation_type[11:13, 15:17] = 1
+
+        nc_water_type[5:9, :5] = 2
+        nc_water_type[11:15, :5] = 2
+
+        nc_soil_type[:, :] = np.where(
+            nc_vegetation_type[:, :] > nc_vegetation_type._FillValue,
+            2,
+            nc_soil_type[:, :])
+        nc_soil_type[:, :] = np.where(
+            nc_pavement_type[:, :] > nc_pavement_type._FillValue,
+            2,
+            nc_soil_type[:, :])
+
+        nc_street_type[:, :] = np.where(nc_pavement_type[:, :] == 1, 11, nc_street_type[:, :])
+        nc_street_type[:, :] = np.where(nc_pavement_type[:, :] == 2, 13, nc_street_type[:, :])
+
+        nc_surface_fraction[0, :, :] = np.where(
+            nc_building_id[:, :] > nc_building_id._FillValue,
+            nc_surface_fraction[0, :, :],
+            0)
+        nc_surface_fraction[2, :, :] = nc_surface_fraction[1, :, :] = nc_surface_fraction[0, :, :]
+        nc_surface_fraction[0, :, :] = np.where(
+            nc_vegetation_type[:, :] > nc_vegetation_type._FillValue,
+            1,
+            nc_surface_fraction[0, :, :])
+        nc_surface_fraction[1, :, :] = np.where(
+            nc_pavement_type[:, :] > nc_pavement_type._FillValue,
+            1,
+            nc_surface_fraction[1, :, :])
+        nc_surface_fraction[2, :, :] = np.where(
+            nc_water_type[:, :] > nc_water_type._FillValue,
+            1,
+            nc_surface_fraction[2, :, :])
+
+        # Set vegetation
+        # --------------
+        lad_profile = [0.0, 0.01070122, 0.1070122, 0.3130108, 0.3879193, 0.1712195]
+        for k in range(len(self.zlad)):
+            nc_lad[k, 7:9, 15:17] = lad_profile[k]
+            nc_lad[k, 11:13, 15:17] = lad_profile[k]
+
     def finalize(self):
-        """ Close file """
+        """Close file."""
         print("Closing file...")
-        
+
         self.nc_file.close()
 
@@ -320,5 +462,4 @@
     driver.write_global_attributes()
     driver.define_dimensions()
-    driver.add_variables()
+    driver.define_variables()
     driver.finalize()
-    
