Changeset 3668 for palm/trunk/SCRIPTS/palm_csd

Timestamp:

Jan 14, 2019 12:49:24 PM (6 years ago)

Author:

maronga

Message:

removed most_methods circular and lookup. added improved version of palm_csd

File:

• palm/trunk/SCRIPTS/palm_csd (modified) (31 diffs)

palm/trunk/SCRIPTS/palm_csd

@@ -25 +25 @@
 # -----------------
 # $Id$
+# Various improvements and bugfixes
+# 
+# 3629 2018-12-13 12:18:54Z maronga
 # Added canopy generator calls. Some improvements
 # 
@@ -72 +75 @@
 
 #  Definition of settings
-   global settings_filename_out
+
    global settings_bridge_width
    global settings_lai_roof_intensive
@@ -100 +103 @@
    global global_site
    global global_source
-   global global_version
+
+   global  path_out
+   global filename_out
+   global version_out
 
 
@@ -158 +164 @@
    global zt_min
 
-   settings_filename_out = "default"
    settings_bridge_width = 3.0
    settings_lai_roof_intensive = 0.5
@@ -185 +190 @@
    global_site = ""
    global_source = ""
-   global_version = 1
-
+
+   path_out = ""
+   version_out = 1
+   filename_out = "default"
+   
    domain_names = []
    domain_px = []
@@ -260 +268 @@
          global_site = config.get(read_tmp, 'site')
          global_source = config.get(read_tmp, 'source')
-         global_version = float(config.get(read_tmp, 'version'))
+
 
       if ( read_tmp == 'settings' ):
-         settings_filename_out = config.get(read_tmp, 'filename_out')
          settings_lai_roof_intensive = config.get(read_tmp, 'lai_roof_intensive')
          settings_lai_roof_extensive = config.get(read_tmp, 'lai_roof_extensive') 
@@ -273 +280 @@
          settings_lai_alpha = float(config.get(read_tmp, 'lai_alpha')) 
          settings_lai_beta = float(config.get(read_tmp, 'lai_beta')) 
+
+      if ( read_tmp == 'output' ):
+         path_out = config.get(read_tmp, 'path')         
+         filename_out = config.get(read_tmp, 'file_out')
+         version_out = float(config.get(read_tmp, 'version'))
 
       if ( read_tmp.split("_")[0] == 'domain' ):
@@ -322 +334 @@
          input_file_street_crossings.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_street_crossings')) 
          input_file_patch_height.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_patch_height')) 
-         input_file_tree_crown_diameter.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_tree_crown_diameter'))  
+         
+         tmp = config.get(read_tmp, 'file_tree_crown_diameter')
+         if tmp is not None:
+            input_file_tree_crown_diameter.append(config.get(read_tmp, 'path') + "/" + tmp)  
+         else:
+            input_file_tree_crown_diameter.append(None)
          input_file_tree_height.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_tree_height'))  
          input_file_tree_trunk_diameter.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_tree_trunk_diameter'))  
@@ -433 +450 @@
 #  Calculate indices and input files
    ii.append(input_px.index(domain_px[i]))
-   filename.append(settings_filename_out + "_" + domain_names[i])
+   filename.append(path_out + "/" + filename_out + "_" + domain_names[i])
    if domain_parent[i] is not None:
       ii_parent.append(domain_names.index(domain_parent[i]))
@@ -447 +464 @@
 #  Create NetCDF output file and set global attributes
    nc_create_file(filename[i])
-   nc_write_global_attributes(filename[i],float(x_UTM[0,0]),float(y_UTM[0,0]),float(lat[0,0]),float(lon[0,0]),"",global_acronym,global_angle,global_author,global_campaign,global_comment,global_contact,global_data_content,global_dependencies,global_institution,global_keywords,global_location,global_palm_version,global_references,global_site,global_source,global_version)
+   nc_write_global_attributes(filename[i],float(x_UTM[0,0]),float(y_UTM[0,0]),float(lat[0,0]),float(lon[0,0]),"",global_acronym,global_angle,global_author,global_campaign,global_comment,global_contact,global_data_content,global_dependencies,global_institution,global_keywords,global_location,global_palm_version,global_references,global_site,global_source,version_out)
 
    
@@ -468 +485 @@
 del zt_all[:]
 
+print( "Shift terrain height by -" + str(zt_min))
 for i in range(0,ndomains):
 
@@ -475 +493 @@
    y = nc_read_from_file_1d(input_file_y[ii[i]], "y", domain_y0[i], domain_y1[i])
 
-   print( "Shift terrain height by -" + str(zt_min))
    zt = zt - zt_min
+   
+   nc_write_global_attribute(filename[i], 'origin_z', float(zt_min))
 
 #  If necessary, interpolate parent domain terrain height on child domain grid and blend the two
@@ -494 +513 @@
 
 #     Interpolate array and bring to PALM grid of child domain
+
       zt_ip = interpolate_2d(zt_parent,tmp_x,tmp_y,x,y)
       zt_ip = bring_to_palm_grid(zt_ip,x,y,domain_dz[ii_parent[i]])
@@ -580 +600 @@
  
    building_type = nc_read_from_file_2d(input_file_building_type[ii[i]], 'Band1', domain_x0[i], domain_x1[i], domain_y0[i], domain_y1[i])  
-   building_type[building_type == 255] = fillvalues["building_type"]
+   building_type[building_type >= 254] = fillvalues["building_type"]
    building_type = np.where(building_type < 1,defaultvalues["building_type"],building_type)
  
@@ -594 +614 @@
       building_type = np.where((building_type == fillvalues["building_type"]) & (buildings_2d != fillvalues["buildings_2d"]),defaultvalues["building_type"],building_type)
       print("Data check #2 " + str(check_arrays_2(buildings_2d,building_type,fillvalues["buildings_2d"],fillvalues["building_type"])))
+ 
  
    nc_write_to_file_2d(filename[i], 'buildings_2d', buildings_2d, datatypes["buildings_2d"],'y','x',fillvalues["buildings_2d"])
@@ -606 +627 @@
    nc_write_attribute(filename[i], 'building_id', 'long_name', 'building id')
    nc_write_attribute(filename[i], 'building_id', 'units', '')
-   nc_write_attribute(filename[i], 'building_id', 'res_orig', domain_px[i]) 
+   nc_write_attribute(filename[i], 'building_id', 'res _orig', domain_px[i]) 
    nc_write_attribute(filename[i], 'building_id', 'coordinates', 'E_UTM N_UTM lon lat') 
    nc_write_attribute(filename[i], 'building_id', 'grid_mapping', 'E_UTM N_UTM lon lat')   
@@ -616 +637 @@
    nc_write_attribute(filename[i], 'building_type', 'coordinates', 'E_UTM N_UTM lon lat') 
    nc_write_attribute(filename[i], 'building_type', 'grid_mapping', 'E_UTM N_UTM lon lat')   
+ 
+   print("out_1: " + str(np.any(building_type == -2)))
  
 del buildings_2d
@@ -705 +728 @@
 #  #6 Remove water for pixels with buildings
    water_type = np.where((water_type != fillvalues["water_type"]) & (building_type != fillvalues["building_type"]),fillvalues["water_type"],water_type) 
-   
-
-#  #7 to be removed: set default soil type everywhere
-   soil_type = np.where((vegetation_type != fillvalues["vegetation_type"]) | (pavement_type != fillvalues["pavement_type"]),defaultvalues["soil_type"],fillvalues["soil_type"]) 
-  
+  
+   
+#  Correct vegetation_type when a vegetation height is available and is indicative of low vegeetation
+   vegetation_height = nc_read_from_file_2d(input_file_vegetation_height[ii[i]], 'Band1', domain_x0[i], domain_x1[i], domain_y0[i], domain_y1[i])   
+         
+   vegetation_type = np.where((vegetation_height != fillvalues["vegetation_height"]) & (vegetation_height == 0.0) & ((vegetation_type == 4) | (vegetation_type == 5) | (vegetation_type == 6) |(vegetation_type == 7) | (vegetation_type == 17) | (vegetation_type == 18)), 3, vegetation_type)
+   vegetation_height = np.where((vegetation_height != fillvalues["vegetation_height"]) & (vegetation_height == 0.0) & ((vegetation_type == 4) | (vegetation_type == 5) | (vegetation_type == 6) |(vegetation_type == 7) | (vegetation_type == 17) | (vegetation_type == 18)), fillvalues["vegetation_height"],vegetation_height)
   
 #  Check for consistency and fill empty fields with default vegetation type
@@ -717 +742 @@
       vegetation_type = np.where(consistency_array == 0,defaultvalues["vegetation_type"],vegetation_type)
       consistency_array, test = check_consistency_4(vegetation_type,building_type,pavement_type,water_type,fillvalues["vegetation_type"],fillvalues["building_type"],fillvalues["pavement_type"],fillvalues["water_type"])
-   
+
+#  #7 to be removed: set default soil type everywhere
+   soil_type = np.where((vegetation_type != fillvalues["vegetation_type"]) | (pavement_type != fillvalues["pavement_type"]),defaultvalues["soil_type"],fillvalues["soil_type"]) 
+        
+   
+#  Check for consistency and fill empty fields with default vegetation type
+   consistency_array, test = check_consistency_3(vegetation_type,pavement_type,soil_type,fillvalues["vegetation_type"],fillvalues["pavement_type"],fillvalues["soil_type"])
+       
 #  Create surface_fraction array
    x = nc_read_from_file_2d_all(filename[i], 'x')
@@ -734 +766 @@
    nc_write_attribute(filename[i], 'surface_fraction', 'res_orig', domain_px[i])     
    del surface_fraction   
-   
- 
-#  Correct vegetation_type when a vegetation height is available and is indicative of low vegeetation
-   vegetation_height = nc_read_from_file_2d(input_file_vegetation_height[ii[i]], 'Band1', domain_x0[i], domain_x1[i], domain_y0[i], domain_y1[i])   
-         
-   vegetation_type = np.where((vegetation_height != fillvalues["vegetation_height"]) & (vegetation_height == 0.0) & ((vegetation_type == 4) | (vegetation_type == 5) | (vegetation_type == 6) |(vegetation_type == 7) | (vegetation_type == 17) | (vegetation_type == 18)), 3, vegetation_type)
-   vegetation_height = np.where((vegetation_height != fillvalues["vegetation_height"]) & (vegetation_height == 0.0) & ((vegetation_type == 4) | (vegetation_type == 5) | (vegetation_type == 6) |(vegetation_type == 7) | (vegetation_type == 17) | (vegetation_type == 18)), fillvalues["vegetation_height"],vegetation_height)
-
 
    nc_write_to_file_2d(filename[i], 'vegetation_type', vegetation_type, datatypes["vegetation_type"],'y','x',fillvalues["vegetation_type"])     
@@ -779 +803 @@
 
 #  pixels with bridges get building_type = 7 = bridge. This does not change the _type setting for the under-bridge
-#  area
+#  area NOTE: when bridges are present the consistency check will fail at the moment
    if domain_3d[i]:
       if np.any(building_type != fillvalues["building_type"]):
@@ -853 +877 @@
 # Read tree data and create LAD and BAD arrays using the canopy generator
 for i in range(0,ndomains):  
+   lai = nc_read_from_file_2d(input_file_lai[ii[i]], 'Band1', domain_x0[i], domain_x1[i], domain_y0[i], domain_y1[i])   
+     
+   vegetation_type = nc_read_from_file_2d_all(filename[i], 'vegetation_type')
+     
+   lai = np.where(vegetation_type == fillvalues["vegetation_type"],fillvalues["vegetation_pars"],lai)
+
+      
+   x = nc_read_from_file_2d_all(filename[i], 'x')
+   y = nc_read_from_file_2d_all(filename[i], 'y')  
+   nvegetation_pars = np.arange(0,12)
+   vegetation_pars = np.ones((len(nvegetation_pars),len(y),len(x)))
+   vegetation_pars[:,:,:] = fillvalues["vegetation_pars"]  
+      
+   vegetation_pars[1,:,:] = lai
+      
+
+#  Write out first version of LAI. Will later be overwritten.
+   nc_write_dimension(filename[i], 'nvegetation_pars', nvegetation_pars, datatypes["nvegetation_pars"]) 
+   nc_write_to_file_3d(filename[i], 'vegetation_pars', vegetation_pars, datatypes["vegetation_pars"],'nvegetation_pars','y','x',fillvalues["vegetation_pars"]) 
+   nc_write_attribute(filename[i], 'vegetation_pars', 'long_name', 'vegetation_pars')
+   nc_write_attribute(filename[i], 'vegetation_pars', 'units', '')
+   nc_write_attribute(filename[i], 'vegetation_pars', 'res_orig', domain_px[i])  
+   nc_write_attribute(filename[i], 'vegetation_pars', 'coordinates', 'E_UTM N_UTM lon lat') 
+   nc_write_attribute(filename[i], 'vegetation_pars', 'grid_mapping', 'E_UTM N_UTM lon lat')   
+
+   del lai, vegetation_pars, vegetation_type
+
+# Read tree data and create LAD and BAD arrays using the canopy generator
+for i in range(0,ndomains):  
    if domain_street_trees[i]: 
-      lai = nc_read_from_file_2d(input_file_lai[ii[i]], 'Band1', domain_x0[i], domain_x1[i], domain_y0[i], domain_y1[i])   
-      
-      
+
+      vegetation_pars =  nc_read_from_file_2d_all(filename[i], 'vegetation_pars')
+
+      lai = np.copy(vegetation_pars[1,:,:])
+
+      x = nc_read_from_file_2d_all(filename[i], 'x')
+      y = nc_read_from_file_2d_all(filename[i], 'y')  
+
 #     Save lai data as default for low and high vegetation
       lai_low = lai
       lai_high = lai
-      
-      x = nc_read_from_file_2d_all(filename[i], 'x')
-      y = nc_read_from_file_2d_all(filename[i], 'y')  
-      nvegetation_pars = np.arange(0,11)
-      vegetation_pars = np.ones((len(nvegetation_pars),len(y),len(x)))
-      vegetation_pars[:,:,:] = fillvalues["vegetation_pars"]  
-      
-      vegetation_pars[1,:,:] = lai
-      
-
-#     Write out first version of LAI. Will later be overwritten.
-      nc_write_dimension(filename[i], 'nvegetation_pars', nvegetation_pars, datatypes["nvegetation_pars"]) 
-      nc_write_to_file_3d(filename[i], 'vegetation_pars', vegetation_pars, datatypes["vegetation_pars"],'nvegetation_pars','y','x',fillvalues["vegetation_pars"]) 
-      nc_write_attribute(filename[i], 'vegetation_pars', 'long_name', 'vegetation_pars')
-      nc_write_attribute(filename[i], 'vegetation_pars', 'units', '')
-      nc_write_attribute(filename[i], 'vegetation_pars', 'res_orig', domain_px[i])  
-      nc_write_attribute(filename[i], 'vegetation_pars', 'coordinates', 'E_UTM N_UTM lon lat') 
-      nc_write_attribute(filename[i], 'vegetation_pars', 'grid_mapping', 'E_UTM N_UTM lon lat')   
-
-
+   
 #     Read all tree parameters from file
       tree_height = nc_read_from_file_2d(input_file_tree_height[ii[i]], 'Band1', domain_x0[i], domain_x1[i], domain_y0[i], domain_y1[i])   
-      tree_crown_diameter = nc_read_from_file_2d(input_file_tree_crown_diameter[ii[i]], 'Band1', domain_x0[i], domain_x1[i], domain_y0[i], domain_y1[i])   
+
+      if (input_file_tree_crown_diameter[ii[i]] is not None):
+         tree_crown_diameter = nc_read_from_file_2d(input_file_tree_crown_diameter[ii[i]], 'Band1', domain_x0[i], domain_x1[i], domain_y0[i], domain_y1[i])   
+         tree_crown_diameter = np.where( (tree_crown_diameter == 0.0) | (tree_crown_diameter == -1.0) ,fillvalues["tree_data"],tree_crown_diameter)
+      else:
+         tree_crown_diameter = np.ones((len(y),len(x)))
+         tree_crown_diameter[:,:] = fillvalues["tree_data"]
+
+
       tree_trunk_diameter = nc_read_from_file_2d(input_file_tree_trunk_diameter[ii[i]], 'Band1', domain_x0[i], domain_x1[i], domain_y0[i], domain_y1[i])         
       tree_type = nc_read_from_file_2d(input_file_tree_type[ii[i]], 'Band1', domain_x0[i], domain_x1[i], domain_y0[i], domain_y1[i]) 
@@ -889 +936 @@
 #     Remove missing values from the data. Reasonable values will be set by the tree generator    
       tree_height = np.where( (tree_height == 0.0) | (tree_height == -1.0) ,fillvalues["tree_data"],tree_height)
-      tree_crown_diameter = np.where( (tree_crown_diameter == 0.0) | (tree_crown_diameter == -1.0) ,fillvalues["tree_data"],tree_crown_diameter)
       tree_trunk_diameter = np.where( (tree_trunk_diameter == 0.0) | (tree_trunk_diameter == -1.0) ,fillvalues["tree_data"],tree_trunk_diameter)
       tree_type = np.where( (tree_type == 0.0) | (tree_type == -1.0) ,fillvalues["tree_data"],tree_type)
@@ -907 +953 @@
 #     For zero-height patches, set vegetation_type to short grass and remove these pixels from the patch height array
       vegetation_type = np.where( (patch_height == 0.0) & ( (vegetation_type == 4) | (vegetation_type == 5) | (vegetation_type == 6) |(vegetation_type == 7) | (vegetation_type == 17) | (vegetation_type == 18) ),3,vegetation_type)
-      patch_type = np.where( (patch_height == 0.0) & ( (vegetation_type == 4) | (vegetation_type == 5) | (vegetation_type == 6) |(vegetation_type == 7) | (vegetation_type == 17) | (vegetation_type == 18) ),fillvalues["tree_data"],patch_height)      
-      nc_overwrite_to_file_2d(filename[i], 'vegetation_type', vegetation_type)  
+      patch_type = np.where( (patch_height == 0.0) & ( (vegetation_type == 4) | (vegetation_type == 5) | (vegetation_type == 6) |(vegetation_type == 7) | (vegetation_type == 17) | (vegetation_type == 18) ),fillvalues["tree_data"],patch_height)   
+      
 
       max_tree_height = max(tree_height.flatten())
@@ -914 +960 @@
       
       if ( (max_tree_height != fillvalues["tree_data"]) | (max_patch_height == fillvalues["tree_data"]) ):
-           
+
          lad, bad, tree_ids, zlad = generate_single_tree_lad(x,y,domain_dz[i],max_tree_height,max_patch_height,tree_type,tree_height,tree_crown_diameter,tree_trunk_diameter,lai,settings_season,fillvalues["tree_data"])
   
@@ -925 +971 @@
             bad[k,:,:] = np.where(buildings_2d == fillvalues["buildings_2d"],bad[k,:,:],fillvalues["tree_data"])
             tree_ids[k,:,:] = np.where(buildings_2d == fillvalues["buildings_2d"],tree_ids[k,:,:],fillvalues["tree_data"])      
-  
+
+         del buildings_2d
   
          nc_write_dimension(filename[i], 'zlad', zlad, datatypes["tree_data"]) 
@@ -949 +996 @@
          nc_write_attribute(filename[i], 'tree_id', 'grid_mapping', 'E_UTM N_UTM lon lat')           
             
-         del buildings_2d, lai, lad, bad, tree_ids, zlad
-
-      del vegetation_pars, tree_height, tree_crown_diameter, tree_trunk_diameter, tree_type, patch_height, vegetation_type, x, y
+         del lai, lad, bad, tree_ids, zlad
+
+      del vegetation_pars, tree_height, tree_crown_diameter, tree_trunk_diameter, tree_type, patch_height, x, y
 
 
@@ -961 +1008 @@
       vegetation_type = nc_read_from_file_2d_all(filename[i], 'vegetation_type') 
       lad = nc_read_from_file_3d_all(filename[i], 'lad') 
+      zlad = nc_read_from_file_1d_all(filename[i], 'zlad') 
       patch_height = nc_read_from_file_2d(input_file_patch_height[ii[i]], 'Band1', domain_x0[i], domain_x1[i], domain_y0[i], domain_y1[i]) 
       vegetation_pars = nc_read_from_file_3d_all(filename[i], 'vegetation_pars') 
       lai = vegetation_pars[1,:,:]
 
-#     For all pixels where the LAD array does not contains a street tree and for which a high vegetation was set and for which a patch height is either missing of larger than one grid spacing
-#     the lai_high is reset to a low vegetation value
-      lai_high = np.where( (lad[0,:,:] == fillvalues["tree_data"]) & ( ( (vegetation_type == 4) | (vegetation_type == 5) | (vegetation_type == 6) |(vegetation_type == 7) | (vegetation_type == 17) | (vegetation_type == 18) ) & ( (patch_height == fillvalues["tree_data"]) | (patch_height >= domain_dz[i])) ),settings_lai_low_default,fillvalues["tree_data"])
-
-#     For all pixels where lai_high is set but no lai is available, set the default high vegetation lai TODO: check if this makes sense bacauce lai_high = lai anway?!
-      lai_high = np.where( (lai_high != fillvalues["tree_data"]) & (lai == fillvalues["tree_data"]), settings_lai_high_default, lai_high)
-
-#     Define a patch height whereever it is missing, but where a high vegetation LAI was set
+
+#     Determine all pixels that do not already have and LAD but which are high vegetation to a dummy value of 1.0 and remove all other pixels
+      lai_high = np.where( (lad[0,:,:] == fillvalues["tree_data"]) & ( ( (vegetation_type == 4) | (vegetation_type == 5) | (vegetation_type == 6) |(vegetation_type == 7) | (vegetation_type == 17) | (vegetation_type == 18) ) & ( (patch_height == fillvalues["tree_data"]) | (patch_height >= domain_dz[i])) ),1.0,fillvalues["tree_data"])
+
+#     Now, assign either the default LAI for high vegetation or the value stored in the LAI array.
+      lai_high = np.where( (lai_high != fillvalues["tree_data"]) & (lai == fillvalues["tree_data"]), settings_lai_high_default, lai)
+
+#     Define a patch height wherever it is missing, but where a high vegetation LAI was set
       patch_height = np.where ( (lai_high != fillvalues["tree_data"]) & (patch_height == fillvalues["tree_data"]), settings_patch_height_default, patch_height)
 
@@ -978 +1026 @@
       patch_height = np.where ( (lad[0,:,:] != fillvalues["tree_data"]), fillvalues["tree_data"], patch_height)
 
+#     Remove patch heights that have no lai_high value
+      patch_height = np.where ( (lai_high == fillvalues["tree_data"]), fillvalues["tree_data"], patch_height)
+
 #     For missing LAI values, set either the high vegetation default or the low vegetation default    
-      lai_high = np.where( (patch_height > 2.0) & (lai_high == fillvalues["tree_data"]),settings_lai_high_default,lai_high)
-      lai_high = np.where( (patch_height <= 2.0) & (lai_high == fillvalues["tree_data"]),settings_lai_low_default,lai_high)
-    
-    
+      lai_high = np.where( (patch_height > 2.0) & (patch_height != fillvalues["tree_data"]) & (lai_high == fillvalues["tree_data"]),settings_lai_high_default,lai_high)
+      lai_high = np.where( (patch_height <= 2.0) & (patch_height != fillvalues["tree_data"]) & (lai_high == fillvalues["tree_data"]),settings_lai_low_default,lai_high)
+         
       if ( max(patch_height.flatten()) >= (2.0 * domain_dz[i]) ):
          print("    start calculating LAD (this might take some time)") 
-         lad_patch, patch_nz, status = process_patch(domain_dz[i],patch_height,lai_high,settings_lai_alpha,settings_lai_beta)
+         lad_patch, patch_nz, status = process_patch(domain_dz[i],patch_height,max(zlad),lai_high,settings_lai_alpha,settings_lai_beta)
 
          lad[0:patch_nz+1,:,:] = np.where( (lad[0:patch_nz+1,:,:] == fillvalues["tree_data"]),lad_patch[0:patch_nz+1,:,:], lad[0:patch_nz+1,:,:] )
@@ -991 +1041 @@
 #     Remove high vegetation wherever it is replaced by a leaf area density. This should effectively remove all high vegetation pixels
       vegetation_type = np.where((lad[0,:,:] != fillvalues["tree_data"]) & (vegetation_type != fillvalues["vegetation_type"]),3,vegetation_type)
-  
- 
+              
 #     If desired, remove all high vegetation. TODO: check if this is still necessary
       if domain_high_vegetation[i]:
@@ -998 +1047 @@
 
   
-#     Remove LAI from pixels where a LAD was set
-      lai_low = np.where( (lad[0,:,:] == fillvalues["tree_data"]), lai,fillvalues["tree_data"])
+#     Set default low LAI for pixels with an LAD (short grass below trees)
+      lai_low = np.where( (lad[0,:,:] == fillvalues["tree_data"]), lai, settings_lai_low_default)
     
-#     Fill low vegetation pixels without LAI set with default value
-      lai_low = np.where( (lai == fillvalues["tree_data"]) & (vegetation_type != fillvalues["vegetation_type"] ), settings_lai_low_default, lai)
-
-    
+#     Fill low vegetation pixels without LAI set or with LAI = 0 with default value
+      lai_low = np.where( ( (lai_low == fillvalues["tree_data"]) | (lai_low == 0.0) ) & (vegetation_type != fillvalues["vegetation_type"] ), settings_lai_low_default, lai_low)
+
+#     Remove lai for pixels that have no vegetation_type
+      lai_low = np.where( vegetation_type != fillvalues["vegetation_type"], lai_low, fillvalues["tree_data"])
+        
 #     Overwrite lai in vegetation_parameters  
-
-      vegetation_pars[1,:,:] = lai_low
+      vegetation_pars[1,:,:] = np.copy(lai_low)
       nc_overwrite_to_file_3d(filename[i], 'vegetation_pars', vegetation_pars) 
 
 #     Overwrite lad array
       nc_overwrite_to_file_3d(filename[i], 'lad', lad) 
-      
-#     Overwrite vegetation_type
-      nc_overwrite_to_file_2d(filename[i], 'vegetation_type', vegetation_type)   
-      
-      
-      del vegetation_type, lad, lai, patch_height, vegetation_pars
+           
+      nc_overwrite_to_file_2d(filename[i], 'vegetation_type', vegetation_type)  
+
+      
+      del vegetation_type, lad, lai, patch_height, vegetation_pars, zlad
+
+# Final consistency check
+for i in range(0,ndomains): 
+   vegetation_type = nc_read_from_file_2d_all(filename[i], 'vegetation_type')  
+   pavement_type = nc_read_from_file_2d_all(filename[i], 'pavement_type')  
+   building_type = nc_read_from_file_2d_all(filename[i], 'building_type')  
+   water_type = nc_read_from_file_2d_all(filename[i], 'water_type')     
+   soil_type = nc_read_from_file_2d_all(filename[i], 'soil_type')  
+
+#  Check for consistency and fill empty fields with default vegetation type
+   consistency_array, test = check_consistency_4(vegetation_type,building_type,pavement_type,water_type,fillvalues["vegetation_type"],fillvalues["building_type"],fillvalues["pavement_type"],fillvalues["water_type"])
+
+#  Check for consistency and fill empty fields with default vegetation type
+   consistency_array, test = check_consistency_3(vegetation_type,pavement_type,soil_type,fillvalues["vegetation_type"],fillvalues["pavement_type"],fillvalues["soil_type"])
+
+   surface_fraction = nc_read_from_file_3d_all(filename[i], 'surface_fraction')       
+   surface_fraction[0,:,:] = np.where(vegetation_type != fillvalues["vegetation_type"], 1.0, 0.0)
+   surface_fraction[1,:,:] = np.where(pavement_type != fillvalues["pavement_type"], 1.0, 0.0)
+   surface_fraction[2,:,:] = np.where(water_type != fillvalues["water_type"], 1.0, 0.0)
+   nc_overwrite_to_file_3d(filename[i], 'surface_fraction', surface_fraction) 
+       
+   del vegetation_type, pavement_type, building_type, water_type, soil_type