Changeset 3668 for palm/trunk

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

Location:

palm/trunk

Files:

• SCRIPTS/.csd.config.default (modified) (10 diffs)
• SCRIPTS/palm_csd (modified) (31 diffs)
• SCRIPTS/palm_csd_files/palm_csd_canopy_generator.py (modified) (11 diffs)
• SCRIPTS/palm_csd_files/palm_csd_netcdf_interface.py (modified) (3 diffs)
• SCRIPTS/palm_csd_files/palm_csd_tools.py (modified) (3 diffs)
• SOURCE/check_parameters.f90 (modified) (2 diffs)
• SOURCE/land_surface_model_mod.f90 (modified) (4 diffs)
• SOURCE/modules.f90 (modified) (2 diffs)
• SOURCE/parin.f90 (modified) (3 diffs)
• SOURCE/read_restart_data_mod.f90 (modified) (3 diffs)
• SOURCE/surface_layer_fluxes_mod.f90 (modified) (9 diffs)
• SOURCE/write_restart_data_mod.f90 (modified) (2 diffs)

palm/trunk/SCRIPTS/.csd.config.default

@@ -25 +25 @@
 # -----------------
 # $Id$
+# Updated configuration
+# 
+# 3567 2018-11-27 13:59:21Z maronga
 # Initial revisions
-#
-#
-# 
-#
 #
 # Description:
@@ -52 +51 @@
 institution = Institute of Meteorology and Climatology, Leibniz University Hannover
 palm_version = 6.0
-version = 1
 rotation_angle = 0.0
 references = text
@@ -58 +56 @@
 
 [settings]
-filename_out = /ldata2/MOSAIK/showcase_python
-lai_season = spring
-lai_high_vegetation_default = 5.0
-lai_low_vegetation_default = 1.0
+lai_roof_intensive = 1.5
+lai_roof_extensive = 3.0
+lai_high_vegetation_default = 6.0
+lai_low_vegetation_default = 3.0
 lai_alpha = 5.0
 lai_beta = 3.0
-tree_height_default = 10.0
+patch_height_default = 10.0
 bridge_width = 3.0
 debug_mode = False
-
+season = "summer"
 
 [input_01]
@@ -79 +77 @@
 file_lon = Berlin_CoordinatesLatLon_x_15m_DLR.nc
 file_zt = Berlin_terrain_height_15m_DLR.nc
-file_buildings_2d = Berlin_building_height_15m_LOD2_DLR.nc
+file_buildings_2d = Berlin_building_height_15m_DLR.nc
 file_building_id = Berlin_building_id_15m_DLR.nc
 file_building_type = Berlin_building_type_15m_DLR.nc
 file_bridges_2d = Berlin_bridges_height_15m_DLR.nc
 file_bridges_id = Berlin_bridges_id_15m_DLR.nc
+file_lai =  Berlin_leaf_area_index_15m_DLR_WANG_summer.nc
 file_vegetation_type = Berlin_vegetation_type_15m_DLR.nc
 file_vegetation_height = Berlin_vegetation_patch_height_15m_DLR.nc
@@ -91 +90 @@
 file_street_type =  Berlin_street_type_15m_DLR.nc
 file_street_crossings =  Berlin_street_crossings_15m_DLR.nc
+file_tree_height =  Berlin_trees_height_clean_15m_DLR.nc
+file_tree_crown_diameter =  Berlin_trees_crown_clean_15m_DLR.nc
+file_tree_trunk_diameter =  Berlin_trees_trunk_clean_15m.nc
+file_tree_type =  Berlin_trees_type_15m_DLR.nc
+file_patch_height =  Berlin_vegetation_patch_height_15m_DLR.nc
+file_vegetation_on_roofs =  Berlin_vegetation_on_roofs_15m_DLR.nc
 
 [input_02]
@@ -102 +107 @@
 file_lon = Berlin_CoordinatesLatLon_x_10m_DLR.nc
 file_zt = Berlin_terrain_height_10m_DLR.nc
-file_buildings_2d = Berlin_building_height_10m_LOD2_DLR.nc
+file_buildings_2d = Berlin_building_height_10m_DLR.nc
 file_building_id = Berlin_building_id_10m_DLR.nc
 file_building_type = Berlin_building_type_10m_DLR.nc
 file_bridges_2d = Berlin_bridges_height_10m_DLR.nc
 file_bridges_id = Berlin_bridges_id_10m_DLR.nc
+file_lai =  Berlin_leaf_area_index_10m_DLR_WANG_summer.nc
 file_vegetation_type = Berlin_vegetation_type_10m_DLR.nc
 file_vegetation_height = Berlin_vegetation_patch_height_10m_DLR.nc
@@ -114 +120 @@
 file_street_type =  Berlin_street_type_10m_DLR.nc
 file_street_crossings =  Berlin_street_crossings_10m_DLR.nc
+file_tree_height =  Berlin_trees_height_clean_10m_DLR.nc
+file_tree_crown_diameter =  Berlin_trees_crown_10m_DLR.nc
+file_tree_trunk_diameter =  Berlin_trees_trunk_clean_10m.nc
+file_tree_type =  Berlin_trees_type_10m_DLR.nc
+file_patch_height =  Berlin_vegetation_patch_height_10m_DLR.nc
+file_vegetation_on_roofs =  Berlin_vegetation_on_roofs_10m_DLR.nc
 
 [input_03]
+path = /ldata2/MOSAIK/Berlin_static_driver_data
+pixel_size = 1.0
 file_x = Berlin_CoordinatesUTM_x_1m_DLR.nc
 file_y = Berlin_CoordinatesUTM_y_1m_DLR.nc
@@ -122 +136 @@
 file_lat = Berlin_CoordinatesLatLon_y_1m_DLR.nc
 file_lon = Berlin_CoordinatesLatLon_x_1m_DLR.nc
-path = /ldata2/MOSAIK/Berlin_static_driver_data
-pixel_size = 1.0
 file_zt = Berlin_terrain_height_1m_DLR.nc
-file_buildings_2d = Berlin_building_height_1m_LOD2_DLR.nc
+file_buildings_2d = Berlin_building_height_1m_DLR.nc
 file_building_id = Berlin_building_id_1m_DLR.nc
 file_bridges_2d = Berlin_bridges_height_1m_DLR.nc
 file_bridges_id = Berlin_bridges_id_1m_DLR.nc
 file_building_type = Berlin_building_type_1m_DLR.nc
+file_lai =  Berlin_leaf_area_index_1m_DLR_WANG_summer.nc
 file_vegetation_type = Berlin_vegetation_type_1m_DLR.nc
 file_vegetation_height = Berlin_vegetation_patch_height_1m_DLR.nc
@@ -137 +150 @@
 file_street_type =  Berlin_street_type_1m_DLR.nc
 file_street_crossings =  Berlin_street_crossings_1m_DLR.nc
+file_tree_height =  Berlin_trees_height_clean_1m_DLR.nc
+file_tree_crown_diameter =  Berlin_trees_crown_clean_1m_DLR.nc
+file_tree_trunk_diameter =  Berlin_trees_trunk_clean_1m.nc
+file_tree_type =  Berlin_trees_type_1m_DLR.nc
+file_patch_height =  Berlin_vegetation_patch_height_1m_DLR.nc
+file_vegetation_on_roofs =  Berlin_vegetation_on_roofs_1m_DLR.nc
+
+[input_04]
+path = /ldata2/MOSAIK/Berlin_static_driver_data
+pixel_size = 2.0
+file_x = Berlin_CoordinatesUTM_x_2m_DLR.nc
+file_y = Berlin_CoordinatesUTM_y_2m_DLR.nc
+file_x_UTM = Berlin_CoordinatesUTM_x_2m_DLR.nc
+file_y_UTM = Berlin_CoordinatesUTM_y_2m_DLR.nc
+file_lat = Berlin_CoordinatesLatLon_y_2m_DLR.nc
+file_lon = Berlin_CoordinatesLatLon_x_2m_DLR.nc
+file_zt = Berlin_terrain_height_2m_DLR.nc
+file_buildings_2d = Berlin_building_height_2m_DLR.nc
+file_building_id = Berlin_building_id_2m_DLR.nc
+file_bridges_2d = Berlin_bridges_height_2m_DLR.nc
+file_bridges_id = Berlin_bridges_id_2m_DLR.nc
+file_building_type = Berlin_building_type_2m_DLR.nc
+file_lai =  Berlin_leaf_area_index_2m_DLR_WANG_summer.nc
+file_vegetation_type = Berlin_vegetation_type_2m_DLR.nc
+file_vegetation_height = Berlin_vegetation_patch_height_2m_DLR.nc
+file_pavement_type = Berlin_pavement_type_2m_DLR.nc
+file_water_type = Berlin_water_type_2m_DLR.nc
+file_soil_type = Berlin_soil_type_2m_DLR.nc
+file_street_type =  Berlin_street_type_2m_DLR.nc
+file_street_crossings =  Berlin_street_crossings_2m_DLR.nc
+file_tree_height =  Berlin_trees_height_clean_2m_DLR.nc
+file_tree_crown_diameter = Berlin_trees_crown_clean_2m_DLR.nc
+file_tree_trunk_diameter =  Berlin_trees_trunk_clean_2m.nc
+file_tree_type =  Berlin_trees_type_2m_DLR.nc
+file_patch_height =  Berlin_vegetation_patch_height_2m_DLR.nc
+file_vegetation_on_roofs =  Berlin_vegetation_on_roofs_2m_DLR.nc
+
 
 [output]
 path = /ldata2/MOSAIK/
-file_out = winter_iop1_test
+file_out = showcase_berlin
 version = 1
 
 [domain_root]
 pixel_size = 15.0
-origin_x = 0.0
-origin_y = 0.0
-nx = 311
-ny = 257
+origin_x = 0
+origin_y = 0
+nx = 3119
+ny = 2573
 buildings_3d = False
 dz = 15.0
-allow_high_vegetation = True
+allow_high_vegetation = False
 generate_vegetation_patches = True
-use_palm_z_axis= True
+use_palm_z_axis= False
 interpolate_terrain = False
 domain_parent
+vegetation_on_roofs = False
+street_trees = True
 
 [domain_N02]
@@ -170 +222 @@
 interpolate_terrain = True
 domain_parent = root
+vegetation_on_roofs = False
+street_trees = True

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

palm/trunk/SCRIPTS/palm_csd_files/palm_csd_canopy_generator.py

@@ -25 +25 @@
 # -----------------
 # $Id$
+# Various improvements and bugfixes
+# 
+# 3629 2018-12-13 12:18:54Z maronga
 # Initial revision
 # 
@@ -72 +75 @@
             if valid_pixels[j,i]:
                tree_id_counter = tree_id_counter + 1
-               print("   Processing tree No " +  str(tree_id_counter) + " ...", end="")
-
+#               print("   Processing tree No " +  str(tree_id_counter) + " ...", end="")
                lad_loc, bad_loc, x_loc, y_loc, z_loc, status = process_single_tree(dx,dz,tree_type[j,i],tree_height[j,i],tree_lai[j,i],tree_dia[j,i],trunk_dia[j,i],season)
                if ( np.any(lad_loc) != fill ):
@@ -101 +103 @@
 
 
-                  if ( status == 0 ):
-                     status_char = " ok."
-                  else:
-                     status_char = " skipped."
-                  print(status_char)
+#                  if ( status == 0 ):
+#                     status_char = " ok."
+#                  else:
+#                     status_char = " skipped."
+#                  print(status_char)
 
                del lad_loc, x_loc, y_loc, z_loc, status
@@ -131 +133 @@
    default_trees = []
                                                 #1   #2   #3   #4    #5   #6   #7   #8     #9
-   default_trees.append(tree("Default",         1.0, 1.0, 4.0, 12.0, 2.0, 0.1, 0.6, 0.025, 0.35)) 
-   default_trees.append(tree("Abies",           3.0, 1.0, 4.0, 12.0, 2.0, 0.1, 0.6, 0.025, 0.80)) 
-   default_trees.append(tree("Acer",            1.0, 1.0, 7.0, 12.0, 2.0, 0.1, 0.6, 0.025, 0.80))
-   default_trees.append(tree("Aesculus",        1.0, 1.0, 7.0, 12.0, 2.0, 0.1, 0.6, 0.025, 1.00))
-   default_trees.append(tree("Ailanthus",       1.0, 1.0, 8.5, 13.5, 2.0, 0.1, 0.6, 0.025, 1.30))
-   default_trees.append(tree("Alnus",           3.0, 1.0, 6.0, 16.0, 2.0, 0.1, 0.6, 0.025, 1.20))
-   default_trees.append(tree("Amelanchier",     1.0, 1.0, 3.0,  4.0, 2.0, 0.1, 0.6, 0.025, 1.20))
-   default_trees.append(tree("Betula",          1.0, 1.0, 6.0, 14.0, 2.0, 0.1, 0.6, 0.025, 0.30))
-   default_trees.append(tree("Buxus",           1.0, 1.0, 4.0,  4.0, 2.0, 0.1, 0.6, 0.025, 0.90))
-   default_trees.append(tree("Calocedrus",      3.0, 1.0, 5.0, 10.0, 2.0, 0.1, 0.6, 0.025, 0.50))
-   default_trees.append(tree("Caragana",        1.0, 1.0, 3.5,  6.0, 2.0, 0.1, 0.6, 0.025, 0.90))
-   default_trees.append(tree("Carpinus",        1.0, 1.0, 6.0, 10.0, 2.0, 0.1, 0.6, 0.025, 0.70))
-   default_trees.append(tree("Carya",           1.0, 1.0, 5.0, 17.0, 2.0, 0.1, 0.6, 0.025, 0.80))
-   default_trees.append(tree("Castanea",        1.0, 1.0, 4.5,  7.0, 2.0, 0.1, 0.6, 0.025, 0.80))
-   default_trees.append(tree("Catalpa",         1.0, 1.0, 5.5,  6.5, 2.0, 0.1, 0.6, 0.025, 0.70))
-   default_trees.append(tree("Cedrus",          1.0, 1.0, 8.0, 13.0, 2.0, 0.1, 0.6, 0.025, 0.80))
-   default_trees.append(tree("Celtis",          1.0, 1.0, 6.0,  9.0, 2.0, 0.1, 0.6, 0.025, 0.80))
-   default_trees.append(tree("Cercidiphyllum",  1.0, 1.0, 3.0,  6.5, 2.0, 0.1, 0.6, 0.025, 0.80))
-   default_trees.append(tree("Cercis",          1.0, 1.0, 2.5,  7.5, 2.0, 0.1, 0.6, 0.025, 0.90))
-   default_trees.append(tree("Chamaecyparis",   5.0, 1.0, 3.5,  9.0, 2.0, 0.1, 0.6, 0.025, 0.70))
-   default_trees.append(tree("Cladrastis",      1.0, 1.0, 5.0, 10.0, 2.0, 0.1, 0.6, 0.025, 0.80))
-   default_trees.append(tree("Cornus",          1.0, 1.0, 4.5,  6.5, 2.0, 0.1, 0.6, 0.025, 1.20))
-   default_trees.append(tree("Corylus",         1.0, 1.0, 5.0,  9.0, 2.0, 0.1, 0.6, 0.025, 0.40))
-   default_trees.append(tree("Cotinus",         1.0, 1.0, 4.0,  4.0, 2.0, 0.1, 0.6, 0.025, 0.70))
-   default_trees.append(tree("Crataegus",       3.0, 1.0, 3.5,  6.0, 2.0, 0.1, 0.6, 0.025, 1.40))
-   default_trees.append(tree("Cryptomeria",     3.0, 1.0, 5.0, 10.0, 2.0, 0.1, 0.6, 0.025, 0.50))
-   default_trees.append(tree("Cupressocyparis", 3.0, 1.0, 3.0,  8.0, 2.0, 0.1, 0.6, 0.025, 0.40))
-   default_trees.append(tree("Cupressus",       3.0, 1.0, 5.0,  7.0, 2.0, 0.1, 0.6, 0.025, 0.40))
-   default_trees.append(tree("Cydonia",         1.0, 1.0, 2.0,  3.0, 2.0, 0.1, 0.6, 0.025, 0.90))
-   default_trees.append(tree("Davidia",         1.0, 1.0,10.0, 14.0, 2.0, 0.1, 0.6, 0.025, 0.40))
-   default_trees.append(tree("Elaeagnus",       1.0, 1.0, 6.5,  6.0, 2.0, 0.1, 0.6, 0.025, 1.20))
-   default_trees.append(tree("Euodia",          1.0, 1.0, 4.5,  6.0, 2.0, 0.1, 0.6, 0.025, 0.90))
-   default_trees.append(tree("Euonymus",        1.0, 1.0, 4.5,  6.0, 2.0, 0.1, 0.6, 0.025, 0.60))
-   default_trees.append(tree("Fagus",           1.0, 1.0,10.0, 12.5, 2.0, 0.1, 0.6, 0.025, 0.50))
-   default_trees.append(tree("Fraxinus",        1.0, 1.0, 5.5, 10.5, 2.0, 0.1, 0.6, 0.025, 1.60))
-   default_trees.append(tree("Ginkgo",          3.0, 1.0, 4.0,  8.5, 2.0, 0.1, 0.6, 0.025, 0.80))
-   default_trees.append(tree("Gleditsia",       1.0, 1.0, 6.5, 10.5, 2.0, 0.1, 0.6, 0.025, 0.60))
-   default_trees.append(tree("Gymnocladus",     1.0, 1.0, 5.5, 10.0, 2.0, 0.1, 0.6, 0.025, 0.80))
-   default_trees.append(tree("Hippophae",       1.0, 1.0, 9.5,  8.5, 2.0, 0.1, 0.6, 0.025, 0.80))
-   default_trees.append(tree("Ilex",            1.0, 1.0, 4.0,  7.5, 2.0, 0.1, 0.6, 0.025, 0.80))
-   default_trees.append(tree("Juglans",         1.0, 1.0, 7.0,  9.0, 2.0, 0.1, 0.6, 0.025, 0.50))
-   default_trees.append(tree("Juniperus",       5.0, 1.0, 3.0,  7.0, 2.0, 0.1, 0.6, 0.025, 0.90))
-   default_trees.append(tree("Koelreuteria",    1.0, 1.0, 3.5,  5.5, 2.0, 0.1, 0.6, 0.025, 0.50))
-   default_trees.append(tree("Laburnum",        1.0, 1.0, 3.0,  6.0, 2.0, 0.1, 0.6, 0.025, 0.60))
-   default_trees.append(tree("Larix",           3.0, 1.0, 7.0, 16.5, 2.0, 0.1, 0.6, 0.025, 0.60))
-   default_trees.append(tree("Ligustrum",       1.0, 1.0, 3.0,  6.0, 2.0, 0.1, 0.6, 0.025, 1.10))
-   default_trees.append(tree("Liquidambar",     3.0, 1.0, 3.0,  7.0, 2.0, 0.1, 0.6, 0.025, 0.30))
-   default_trees.append(tree("Liriodendron",    3.0, 1.0, 4.5,  9.5, 2.0, 0.1, 0.6, 0.025, 0.50))
-   default_trees.append(tree("Lonicera",        1.0, 1.0, 7.0,  9.0, 2.0, 0.1, 0.6, 0.025, 0.70))
-   default_trees.append(tree("Magnolia",        1.0, 1.0, 3.0,  5.0, 2.0, 0.1, 0.6, 0.025, 0.60))
-   default_trees.append(tree("Malus",           1.0, 1.0, 4.5,  5.0, 2.0, 0.1, 0.6, 0.025, 0.30))
-   default_trees.append(tree("Metasequoia",     5.0, 1.0, 4.5, 12.0, 2.0, 0.1, 0.6, 0.025, 0.50))
-   default_trees.append(tree("Morus",           1.0, 1.0, 7.5, 11.5, 2.0, 0.1, 0.6, 0.025, 1.00))
-   default_trees.append(tree("Ostrya",          1.0, 1.0, 2.0,  6.0, 2.0, 0.1, 0.6, 0.025, 1.00))
-   default_trees.append(tree("Parrotia",        1.0, 1.0, 7.0,  7.0, 2.0, 0.1, 0.6, 0.025, 0.30))
-   default_trees.append(tree("Paulownia",       1.0, 1.0, 4.0,  8.0, 2.0, 0.1, 0.6, 0.025, 0.40))
-   default_trees.append(tree("Phellodendron",   1.0, 1.0,13.5, 13.5, 2.0, 0.1, 0.6, 0.025, 0.50))
-   default_trees.append(tree("Picea",           3.0, 1.0, 3.0, 13.0, 2.0, 0.1, 0.6, 0.025, 0.90))
-   default_trees.append(tree("Pinus",           3.0, 1.0, 6.0, 16.0, 2.0, 0.1, 0.6, 0.025, 0.80))
-   default_trees.append(tree("Platanus",        1.0, 1.0,10.0, 14.5, 2.0, 0.1, 0.6, 0.025, 1.10))
-   default_trees.append(tree("Populus",         1.0, 1.0, 9.0, 20.0, 2.0, 0.1, 0.6, 0.025, 1.40))
-   default_trees.append(tree("Prunus",          1.0, 1.0, 5.0,  7.0, 2.0, 0.1, 0.6, 0.025, 1.60))
-   default_trees.append(tree("Pseudotsuga",     3.0, 1.0, 6.0, 17.5, 2.0, 0.1, 0.6, 0.025, 0.70))
-   default_trees.append(tree("Ptelea",          1.0, 1.0, 5.0,  4.0, 2.0, 0.1, 0.6, 0.025, 1.10))
-   default_trees.append(tree("Pterocaria",      1.0, 1.0,10.0, 12.0, 2.0, 0.1, 0.6, 0.025, 0.50))
-   default_trees.append(tree("Pterocarya",      1.0, 1.0,11.5, 14.5, 2.0, 0.1, 0.6, 0.025, 1.60))
-   default_trees.append(tree("Pyrus",           3.0, 1.0, 3.0,  6.0, 2.0, 0.1, 0.6, 0.025, 1.80))
+   default_trees.append(tree("Default",         1.0, 1.0, 4.0, 12.0, 3.0, 0.8, 0.6, 0.025, 0.35)) 
+   default_trees.append(tree("Abies",           3.0, 1.0, 4.0, 12.0, 3.0, 0.8, 0.6, 0.025, 0.80)) 
+   default_trees.append(tree("Acer",            1.0, 1.0, 7.0, 12.0, 3.0, 0.8, 0.6, 0.025, 0.80))
+   default_trees.append(tree("Aesculus",        1.0, 1.0, 7.0, 12.0, 3.0, 0.8, 0.6, 0.025, 1.00))
+   default_trees.append(tree("Ailanthus",       1.0, 1.0, 8.5, 13.5, 3.0, 0.8, 0.6, 0.025, 1.30))
+   default_trees.append(tree("Alnus",           3.0, 1.0, 6.0, 16.0, 3.0, 0.8, 0.6, 0.025, 1.20))
+   default_trees.append(tree("Amelanchier",     1.0, 1.0, 3.0,  4.0, 3.0, 0.8, 0.6, 0.025, 1.20))
+   default_trees.append(tree("Betula",          1.0, 1.0, 6.0, 14.0, 3.0, 0.8, 0.6, 0.025, 0.30))
+   default_trees.append(tree("Buxus",           1.0, 1.0, 4.0,  4.0, 3.0, 0.8, 0.6, 0.025, 0.90))
+   default_trees.append(tree("Calocedrus",      3.0, 1.0, 5.0, 10.0, 3.0, 0.8, 0.6, 0.025, 0.50))
+   default_trees.append(tree("Caragana",        1.0, 1.0, 3.5,  6.0, 3.0, 0.8, 0.6, 0.025, 0.90))
+   default_trees.append(tree("Carpinus",        1.0, 1.0, 6.0, 10.0, 3.0, 0.8, 0.6, 0.025, 0.70))
+   default_trees.append(tree("Carya",           1.0, 1.0, 5.0, 17.0, 3.0, 0.8, 0.6, 0.025, 0.80))
+   default_trees.append(tree("Castanea",        1.0, 1.0, 4.5,  7.0, 3.0, 0.8, 0.6, 0.025, 0.80))
+   default_trees.append(tree("Catalpa",         1.0, 1.0, 5.5,  6.5, 3.0, 0.8, 0.6, 0.025, 0.70))
+   default_trees.append(tree("Cedrus",          1.0, 1.0, 8.0, 13.0, 3.0, 0.8, 0.6, 0.025, 0.80))
+   default_trees.append(tree("Celtis",          1.0, 1.0, 6.0,  9.0, 3.0, 0.8, 0.6, 0.025, 0.80))
+   default_trees.append(tree("Cercidiphyllum",  1.0, 1.0, 3.0,  6.5, 3.0, 0.8, 0.6, 0.025, 0.80))
+   default_trees.append(tree("Cercis",          1.0, 1.0, 2.5,  7.5, 3.0, 0.8, 0.6, 0.025, 0.90))
+   default_trees.append(tree("Chamaecyparis",   5.0, 1.0, 3.5,  9.0, 3.0, 0.8, 0.6, 0.025, 0.70))
+   default_trees.append(tree("Cladrastis",      1.0, 1.0, 5.0, 10.0, 3.0, 0.8, 0.6, 0.025, 0.80))
+   default_trees.append(tree("Cornus",          1.0, 1.0, 4.5,  6.5, 3.0, 0.8, 0.6, 0.025, 1.20))
+   default_trees.append(tree("Corylus",         1.0, 1.0, 5.0,  9.0, 3.0, 0.8, 0.6, 0.025, 0.40))
+   default_trees.append(tree("Cotinus",         1.0, 1.0, 4.0,  4.0, 3.0, 0.8, 0.6, 0.025, 0.70))
+   default_trees.append(tree("Crataegus",       3.0, 1.0, 3.5,  6.0, 3.0, 0.8, 0.6, 0.025, 1.40))
+   default_trees.append(tree("Cryptomeria",     3.0, 1.0, 5.0, 10.0, 3.0, 0.8, 0.6, 0.025, 0.50))
+   default_trees.append(tree("Cupressocyparis", 3.0, 1.0, 3.0,  8.0, 3.0, 0.8, 0.6, 0.025, 0.40))
+   default_trees.append(tree("Cupressus",       3.0, 1.0, 5.0,  7.0, 3.0, 0.8, 0.6, 0.025, 0.40))
+   default_trees.append(tree("Cydonia",         1.0, 1.0, 2.0,  3.0, 3.0, 0.8, 0.6, 0.025, 0.90))
+   default_trees.append(tree("Davidia",         1.0, 1.0,10.0, 14.0, 3.0, 0.8, 0.6, 0.025, 0.40))
+   default_trees.append(tree("Elaeagnus",       1.0, 1.0, 6.5,  6.0, 3.0, 0.8, 0.6, 0.025, 1.20))
+   default_trees.append(tree("Euodia",          1.0, 1.0, 4.5,  6.0, 3.0, 0.8, 0.6, 0.025, 0.90))
+   default_trees.append(tree("Euonymus",        1.0, 1.0, 4.5,  6.0, 3.0, 0.8, 0.6, 0.025, 0.60))
+   default_trees.append(tree("Fagus",           1.0, 1.0,10.0, 12.5, 3.0, 0.8, 0.6, 0.025, 0.50))
+   default_trees.append(tree("Fraxinus",        1.0, 1.0, 5.5, 10.5, 3.0, 0.8, 0.6, 0.025, 1.60))
+   default_trees.append(tree("Ginkgo",          3.0, 1.0, 4.0,  8.5, 3.0, 0.8, 0.6, 0.025, 0.80))
+   default_trees.append(tree("Gleditsia",       1.0, 1.0, 6.5, 10.5, 3.0, 0.8, 0.6, 0.025, 0.60))
+   default_trees.append(tree("Gymnocladus",     1.0, 1.0, 5.5, 10.0, 3.0, 0.8, 0.6, 0.025, 0.80))
+   default_trees.append(tree("Hippophae",       1.0, 1.0, 9.5,  8.5, 3.0, 0.8, 0.6, 0.025, 0.80))
+   default_trees.append(tree("Ilex",            1.0, 1.0, 4.0,  7.5, 3.0, 0.8, 0.6, 0.025, 0.80))
+   default_trees.append(tree("Juglans",         1.0, 1.0, 7.0,  9.0, 3.0, 0.8, 0.6, 0.025, 0.50))
+   default_trees.append(tree("Juniperus",       5.0, 1.0, 3.0,  7.0, 3.0, 0.8, 0.6, 0.025, 0.90))
+   default_trees.append(tree("Koelreuteria",    1.0, 1.0, 3.5,  5.5, 3.0, 0.8, 0.6, 0.025, 0.50))
+   default_trees.append(tree("Laburnum",        1.0, 1.0, 3.0,  6.0, 3.0, 0.8, 0.6, 0.025, 0.60))
+   default_trees.append(tree("Larix",           3.0, 1.0, 7.0, 16.5, 3.0, 0.8, 0.6, 0.025, 0.60))
+   default_trees.append(tree("Ligustrum",       1.0, 1.0, 3.0,  6.0, 3.0, 0.8, 0.6, 0.025, 1.10))
+   default_trees.append(tree("Liquidambar",     3.0, 1.0, 3.0,  7.0, 3.0, 0.8, 0.6, 0.025, 0.30))
+   default_trees.append(tree("Liriodendron",    3.0, 1.0, 4.5,  9.5, 3.0, 0.8, 0.6, 0.025, 0.50))
+   default_trees.append(tree("Lonicera",        1.0, 1.0, 7.0,  9.0, 3.0, 0.8, 0.6, 0.025, 0.70))
+   default_trees.append(tree("Magnolia",        1.0, 1.0, 3.0,  5.0, 3.0, 0.8, 0.6, 0.025, 0.60))
+   default_trees.append(tree("Malus",           1.0, 1.0, 4.5,  5.0, 3.0, 0.8, 0.6, 0.025, 0.30))
+   default_trees.append(tree("Metasequoia",     5.0, 1.0, 4.5, 12.0, 3.0, 0.8, 0.6, 0.025, 0.50))
+   default_trees.append(tree("Morus",           1.0, 1.0, 7.5, 11.5, 3.0, 0.8, 0.6, 0.025, 1.00))
+   default_trees.append(tree("Ostrya",          1.0, 1.0, 2.0,  6.0, 3.0, 0.8, 0.6, 0.025, 1.00))
+   default_trees.append(tree("Parrotia",        1.0, 1.0, 7.0,  7.0, 3.0, 0.8, 0.6, 0.025, 0.30))
+   default_trees.append(tree("Paulownia",       1.0, 1.0, 4.0,  8.0, 3.0, 0.8, 0.6, 0.025, 0.40))
+   default_trees.append(tree("Phellodendron",   1.0, 1.0,13.5, 13.5, 3.0, 0.8, 0.6, 0.025, 0.50))
+   default_trees.append(tree("Picea",           3.0, 1.0, 3.0, 13.0, 3.0, 0.8, 0.6, 0.025, 0.90))
+   default_trees.append(tree("Pinus",           3.0, 1.0, 6.0, 16.0, 3.0, 0.8, 0.6, 0.025, 0.80))
+   default_trees.append(tree("Platanus",        1.0, 1.0,10.0, 14.5, 3.0, 0.8, 0.6, 0.025, 1.10))
+   default_trees.append(tree("Populus",         1.0, 1.0, 9.0, 20.0, 3.0, 0.8, 0.6, 0.025, 1.40))
+   default_trees.append(tree("Prunus",          1.0, 1.0, 5.0,  7.0, 3.0, 0.8, 0.6, 0.025, 1.60))
+   default_trees.append(tree("Pseudotsuga",     3.0, 1.0, 6.0, 17.5, 3.0, 0.8, 0.6, 0.025, 0.70))
+   default_trees.append(tree("Ptelea",          1.0, 1.0, 5.0,  4.0, 3.0, 0.8, 0.6, 0.025, 1.10))
+   default_trees.append(tree("Pterocaria",      1.0, 1.0,10.0, 12.0, 3.0, 0.8, 0.6, 0.025, 0.50))
+   default_trees.append(tree("Pterocarya",      1.0, 1.0,11.5, 14.5, 3.0, 0.8, 0.6, 0.025, 1.60))
+   default_trees.append(tree("Pyrus",           3.0, 1.0, 3.0,  6.0, 3.0, 0.8, 0.6, 0.025, 1.80))
    default_trees.append(tree("Quercus",         1.0, 1.0, 8.0, 14.0, 3.1, 0.1, 0.6, 0.025, 0.40)) #
-   default_trees.append(tree("Rhamnus",         1.0, 1.0, 4.5,  4.5, 2.0, 0.1, 0.6, 0.025, 1.30))
-   default_trees.append(tree("Rhus",            1.0, 1.0, 7.0,  5.5, 2.0, 0.1, 0.6, 0.025, 0.50))
-   default_trees.append(tree("Robinia",         1.0, 1.0, 4.5, 13.5, 2.0, 0.1, 0.6, 0.025, 0.50))
-   default_trees.append(tree("Salix",           1.0, 1.0, 7.0, 14.0, 2.0, 0.1, 0.6, 0.025, 1.10))
-   default_trees.append(tree("Sambucus",        1.0, 1.0, 8.0,  6.0, 2.0, 0.1, 0.6, 0.025, 1.40))
-   default_trees.append(tree("Sasa",            1.0, 1.0,10.0, 25.0, 2.0, 0.1, 0.6, 0.025, 0.60))
-   default_trees.append(tree("Sequoiadendron",  5.0, 1.0, 5.5, 10.5, 2.0, 0.1, 0.6, 0.025, 1.60))
-   default_trees.append(tree("Sophora",         1.0, 1.0, 7.5, 10.0, 2.0, 0.1, 0.6, 0.025, 1.40))
-   default_trees.append(tree("Sorbus",          1.0, 1.0, 4.0,  7.0, 2.0, 0.1, 0.6, 0.025, 1.10))
-   default_trees.append(tree("Syringa",         1.0, 1.0, 4.5,  5.0, 2.0, 0.1, 0.6, 0.025, 0.60))
-   default_trees.append(tree("Tamarix",         1.0, 1.0, 6.0,  7.0, 2.0, 0.1, 0.6, 0.025, 0.50))
-   default_trees.append(tree("Taxodium",        5.0, 1.0, 6.0, 16.5, 2.0, 0.1, 0.6, 0.025, 0.60))
-   default_trees.append(tree("Taxus",           2.0, 1.0, 5.0,  7.5, 2.0, 0.1, 0.6, 0.025, 1.50))
-   default_trees.append(tree("Thuja",           3.0, 1.0, 3.5,  9.0, 2.0, 0.1, 0.6, 0.025, 0.70))
-   default_trees.append(tree("Tilia",           3.0, 1.0, 7.0, 12.5, 2.0, 0.1, 0.6, 0.025, 0.70))
-   default_trees.append(tree("Tsuga",           3.0, 1.0, 6.0, 10.5, 2.0, 0.1, 0.6, 0.025, 1.10))
-   default_trees.append(tree("Ulmus",           1.0, 1.0, 7.5, 14.0, 2.0, 0.1, 0.6, 0.025, 0.80))
-   default_trees.append(tree("Zelkova",         1.0, 1.0, 4.0,  5.5, 2.0, 0.1, 0.6, 0.025, 1.20))
-   default_trees.append(tree("Zenobia",         1.0, 1.0, 5.0,  5.0, 2.0, 0.1, 0.6, 0.025, 0.40))
+   default_trees.append(tree("Rhamnus",         1.0, 1.0, 4.5,  4.5, 3.0, 0.8, 0.6, 0.025, 1.30))
+   default_trees.append(tree("Rhus",            1.0, 1.0, 7.0,  5.5, 3.0, 0.8, 0.6, 0.025, 0.50))
+   default_trees.append(tree("Robinia",         1.0, 1.0, 4.5, 13.5, 3.0, 0.8, 0.6, 0.025, 0.50))
+   default_trees.append(tree("Salix",           1.0, 1.0, 7.0, 14.0, 3.0, 0.8, 0.6, 0.025, 1.10))
+   default_trees.append(tree("Sambucus",        1.0, 1.0, 8.0,  6.0, 3.0, 0.8, 0.6, 0.025, 1.40))
+   default_trees.append(tree("Sasa",            1.0, 1.0,10.0, 25.0, 3.0, 0.8, 0.6, 0.025, 0.60))
+   default_trees.append(tree("Sequoiadendron",  5.0, 1.0, 5.5, 10.5, 3.0, 0.8, 0.6, 0.025, 1.60))
+   default_trees.append(tree("Sophora",         1.0, 1.0, 7.5, 10.0, 3.0, 0.8, 0.6, 0.025, 1.40))
+   default_trees.append(tree("Sorbus",          1.0, 1.0, 4.0,  7.0, 3.0, 0.8, 0.6, 0.025, 1.10))
+   default_trees.append(tree("Syringa",         1.0, 1.0, 4.5,  5.0, 3.0, 0.8, 0.6, 0.025, 0.60))
+   default_trees.append(tree("Tamarix",         1.0, 1.0, 6.0,  7.0, 3.0, 0.8, 0.6, 0.025, 0.50))
+   default_trees.append(tree("Taxodium",        5.0, 1.0, 6.0, 16.5, 3.0, 0.8, 0.6, 0.025, 0.60))
+   default_trees.append(tree("Taxus",           2.0, 1.0, 5.0,  7.5, 3.0, 0.8, 0.6, 0.025, 1.50))
+   default_trees.append(tree("Thuja",           3.0, 1.0, 3.5,  9.0, 3.0, 0.8, 0.6, 0.025, 0.70))
+   default_trees.append(tree("Tilia",           3.0, 1.0, 7.0, 12.5, 3.0, 0.8, 0.6, 0.025, 0.70))
+   default_trees.append(tree("Tsuga",           3.0, 1.0, 6.0, 10.5, 3.0, 0.8, 0.6, 0.025, 1.10))
+   default_trees.append(tree("Ulmus",           1.0, 1.0, 7.5, 14.0, 3.0, 0.8, 0.6, 0.025, 0.80))
+   default_trees.append(tree("Zelkova",         1.0, 1.0, 4.0,  5.5, 3.0, 0.8, 0.6, 0.025, 1.20))
+   default_trees.append(tree("Zenobia",         1.0, 1.0, 5.0,  5.0, 3.0, 0.8, 0.6, 0.025, 0.40))
    
 #  Define fill values
@@ -237 +239 @@
       if (season == "summer"):
          tree_lai = default_trees[tree_type].lai_summer
+
       else:
          tree_lai = default_trees[tree_type].lai_winter      
@@ -279 +282 @@
    
    lad_max = tree_lai / (lad_max_part_1[0] + lad_max_part_2[0])
-
    
    
@@ -298 +300 @@
    z = np.arange(0,nz*dz,dz)
    z[1:] = z[1:] - 0.5 * dz
-
 
 #  Define center of the tree position inside the local LAD domain
@@ -468 +469 @@
 
 
-def process_patch(dz,patch_height,patch_lai,alpha,beta):
+def process_patch(dz,patch_height,max_height_lad,patch_lai,alpha,beta):
    
 #  Define fill values
@@ -481 +482 @@
    pch_index = np.where( (pch_index[:,:] == -1) ,0,pch_index[:,:])
 
-   max_canopy_height = max(patch_height.flatten())
+   max_canopy_height = max(max(patch_height.flatten()),max_height_lad)
 
    z = np.arange(0,math.floor(max_canopy_height/dz)*dz+2*dz,dz)
@@ -500 +501 @@
           int_bpdf = 0.0
           if ( (patch_height[j,i] != fillvalues["tree_data"]) & (patch_height[j,i] >= (0.5*dz)) ):
-             for k in range(0,pch_index[j,i]):
+             for k in range(1,pch_index[j,i]):
                 int_bpdf = int_bpdf + ( ( ( z[k] / patch_height[j,i] )**( alpha - 1 ) ) * ( ( 1.0 - ( z[k] / patch_height[j,i] ) )**(beta - 1 ) ) * ( dz / patch_height[j,i] ) )
 
-             for k in range(0,pch_index[j,i]):
+             for k in range(1,pch_index[j,i]):
                 pre_lad[k] =  patch_lai[j,i] * ( ( ( dz*k / patch_height[j,i] )**( alpha - 1.0 ) ) * ( ( 1.0 - ( dz*k / patch_height[j,i] ) )**(beta - 1.0 ) ) / int_bpdf ) / patch_height[j,i]
 
@@ -538 +539 @@
 # dbh: default trunk diameter at breast height (1.4 m) (m)
 #
-#
 class tree:
     def __init__(self, species, shape, ratio, diameter, height, lai_summer, lai_winter, lad_max_height, bad_scale, dbh): 

palm/trunk/SCRIPTS/palm_csd_files/palm_csd_netcdf_interface.py

@@ -25 +25 @@
 # -----------------
 # $Id$
+# Some improvements and new routines
+# 
+# 3629 2018-12-13 12:18:54Z maronga
 # Some new routines
 # 
@@ -39 +42 @@
 
 from netCDF4 import Dataset
+
+def nc_read_from_file_1d_all(filename, varname):
+
+
+   import numpy as np
+   import sys
+   
+   try:
+      f = open(filename)
+      f.close()
+#      print("Load: " + filename + ".")
+   except FileNotFoundError:
+      print("Error: " + filename + ". No such file. Aborting...")
+      sys.exit(1)
+   
+   nc_file = Dataset(filename, "r+", format="NETCDF4")
+   tmp_array = np.array(nc_file.variables[varname][:] , dtype=type(nc_file.variables[varname]))
+
+   return tmp_array
 
 def nc_read_from_file_2d(filename, varname, x0, x1, y0, y1):
@@ -174 +196 @@
    return 0
 
+def nc_write_global_attribute(filename,attribute,value):
+
+   import datetime
+   import os
+
+   print("Writing attribute " + attribute + " to file...")
+
+   f =  Dataset(filename, "a", format="NETCDF4")
+   
+   f.setncattr(attribute,value)
+
+   f.close()
+
 
 def nc_write_dimension(filename,varname,array,datatype):

palm/trunk/SCRIPTS/palm_csd_files/palm_csd_tools.py

@@ -25 +25 @@
 # -----------------
 # $Id$
+# Minor changes
+# 
+# 3567 2018-11-27 13:59:21Z maronga
 # Initial revisions
-#
-#
-# 
 #
 #
@@ -77 +77 @@
 
 def interpolate_2d(array,x1,y1,x2,y2):
-    
-      tmp_int2d = interp2d(y1,x1,array,kind='linear')
-      array_ip = tmp_int2d(y2.astype(float), x2.astype(float)) 
+   tmp_int2d = interp2d(x1,y1,array,kind='linear')
+   array_ip = tmp_int2d(x2.astype(float), y2.astype(float)) 
 
-      return array_ip 
+   return array_ip 
    
    
@@ -136 +135 @@
    return result
 
+def check_consistency_3(array1,array2,array3,fill1,fill2,fill3):
+
+   tmp_array = np.where(array1 != fill1,1,0) + np.where(array2 != fill2,1,0) + np.where(array3 != fill3,-1,0)
+   
+   test = np.any(tmp_array.flatten() != 0)
+   if test:
+      print("soil_type array is not consistent!")
+      print("max: " + str(max(tmp_array.flatten())) + ", min: " + str(min(tmp_array.flatten())))
+   else:
+      print("soil_type array is consistent!")     
+   return tmp_array, test
+
+
+
 
 def check_consistency_4(array1,array2,array3,array4,fill1,fill2,fill3,fill4):

palm/trunk/SOURCE/check_parameters.f90

@@ -25 +25 @@
 ! -----------------
 ! $Id$
+! most_method removed
+! 
+! 3655 2019-01-07 16:51:22Z knoop
 ! Formatting
 ! 
@@ -3790 +3793 @@
     ENDIF
 
-!
-!-- Check for valid setting of most_method
-    IF ( TRIM( most_method ) /= 'circular'  .AND.                              &
-         TRIM( most_method ) /= 'newton'    .AND.                              &
-         TRIM( most_method ) /= 'lookup' )  THEN
-       message_string = 'most_method = "' // TRIM( most_method ) //            &
-                        '" is unknown'
-       CALL message( 'check_parameters', 'PA0416', 1, 2, 0, 6, 0 )
-    ENDIF
 
 !

palm/trunk/SOURCE/land_surface_model_mod.f90

@@ -25 +25 @@
 ! -----------------
 ! $Id$
+! Removed most_method
+! 
+! 3655 2019-01-07 16:51:22Z knoop
 ! nopointer option removed
 ! 
@@ -1324 +1327 @@
 
     USE control_parameters,                                                    &
-        ONLY:  bc_pt_b, bc_q_b, constant_flux_layer, message_string,           &
-               most_method
+        ONLY:  bc_pt_b, bc_q_b, constant_flux_layer, message_string
                      
     
@@ -1446 +1448 @@
     
     IF ( TRIM( surface_type ) == 'water' )  THEN
-
-       IF ( TRIM( most_method ) == 'lookup' )  THEN    
-          WRITE( message_string, * ) 'surface_type = ', surface_type,          &
-                                     ' is not allowed in combination with ',   &
-                                     'most_method = ', most_method
-          CALL message( 'lsm_check_parameters', 'PA0414', 1, 2, 0, 6, 0 )
-       ENDIF
 
        IF ( water_type == 0 )  THEN  
@@ -1531 +1526 @@
 
     IF ( TRIM( surface_type ) == 'netcdf' )  THEN
-       IF ( ANY( water_type_f%var /= water_type_f%fill )  .AND.                &
-            TRIM( most_method ) == 'lookup' )  THEN    
-          WRITE( message_string, * ) 'water-surfaces are not allowed in ' //   &
-                                     'combination with most_method = ',        &
-                                     TRIM( most_method )
-          CALL message( 'lsm_check_parameters', 'PA0999', 2, 2, 0, 6, 0 )
-       ENDIF
 !
 !--    MS: Some problme here, after calling message everythings stucks at

palm/trunk/SOURCE/modules.f90

@@ -25 +25 @@
 ! -----------------
 ! $Id$
+! Removed most_method
+! 
+! 3648 2019-01-02 16:35:46Z suehring
 ! -surface_data_output +surface_output
 ! 
@@ -1090 +1093 @@
     CHARACTER (LEN=8)    ::  coupling_char = ''                           !< appended to filenames in coupled or nested runs ('_O': ocean PE,
                                                                           !< '_NV': vertically nested atmosphere PE, '_N##': PE of nested domain ##
-    CHARACTER (LEN=8)    ::  most_method = 'newton'                       !< namelist parameter
     CHARACTER (LEN=10)   ::  run_date = ' '                               !< date of simulation run
     CHARACTER (LEN=8)    ::  run_time = ' '                               !< time of simulation run

palm/trunk/SOURCE/parin.f90

@@ -25 +25 @@
 ! -----------------
 ! $Id$
+! Removed most_method
+! 
+! 3649 2019-01-02 16:52:21Z suehring
 ! Delete debug-print statements
 ! 
@@ -573 +576 @@
              loop_optimization, lsf_exception, masking_method, mg_cycles,      &
              mg_switch_to_pe0_level, mixing_length_1d, momentum_advec,         &
-             most_method,                                                      &
              netcdf_precision, neutral, ngsrb,                                 &
              nsor, nsor_ini, nudging, nx, ny, nz, ocean_mode, omega,           &
@@ -644 +646 @@
              loop_optimization, lsf_exception, masking_method, mg_cycles,      &
              mg_switch_to_pe0_level, mixing_length_1d, momentum_advec,         &
-             most_method,                                                      &
              netcdf_precision, neutral, ngsrb,                                 &
              nsor, nsor_ini, nudging, nx, ny, nz, ocean_mode, omega,           &

palm/trunk/SOURCE/read_restart_data_mod.f90

@@ -25 +25 @@
 ! -----------------
 ! $Id$
+! Removed most_method and increased binary version
+! 
+! 3655 2019-01-07 16:51:22Z knoop
 ! Implementation of the PALM module interface
 ! 
@@ -203 +206 @@
        READ ( 13 )  version_on_file
 
-       binary_version_global = '4.7'
+       binary_version_global = '4.8'
        IF ( TRIM( version_on_file ) /= TRIM( binary_version_global ) )  THEN
           WRITE( message_string, * ) 'version mismatch concerning ',           &
@@ -488 +491 @@
              CASE ( 'momentum_advec' )
                 READ ( 13 )  momentum_advec
-             CASE ( 'most_method' )
-                READ ( 13 )  most_method
              CASE ( 'netcdf_precision' )
                 READ ( 13 )  netcdf_precision

palm/trunk/SOURCE/surface_layer_fluxes_mod.f90

@@ -26 +26 @@
 ! -----------------
 ! $Id$
+! Removed methods "circular" and "lookup"
+! 
+! 3655 2019-01-07 16:51:22Z knoop
 ! OpenACC port for SPEC
 ! 
@@ -228 +231 @@
 ! ------------
 !> Diagnostic computation of vertical fluxes in the constant flux layer from the
-!> values of the variables at grid point k=1. Three different methods are
-!> available:
-!> 1) the "old" version (most_method = 'circular') which is fast, but inaccurate
-!> 2) a Newton iteration method (most_method = 'newton'), which is accurate, but
-!>    slower
-!> 3) a method using a lookup table which is fast and accurate. Note, however,
-!>    that this method cannot be used in case of roughness heterogeneity
+!> values of the variables at grid point k=1 based on Newton iteration
 !>
 !> @todo (re)move large_scale_forcing actions
@@ -261 +258 @@
                intermediate_timestep_count, intermediate_timestep_count_max,   &
                land_surface, large_scale_forcing, lsf_surf, message_string,    &
-               most_method, neutral, passive_scalar, pt_surface, q_surface,    &
+               neutral, passive_scalar, pt_surface, q_surface,                 &
                run_coupled, surface_pressure, simulated_time, terminate_run,   &
                time_since_reference_point, urban_surface,                      &
@@ -400 +397 @@
 !--    First, calculate the new Obukhov length, then new friction velocity,
 !--    followed by the new scaling parameters (th*, q*, etc.), and the new
-!--    surface fluxes if required. The old routine ("circular") requires a 
-!--    different order of calls as the scaling parameters from the previous time
-!--    steps are used to calculate the Obukhov length
-
-!
-!--    Depending on setting of most_method use the "old" routine
-!--    Note, each routine is called for different surface types. 
+!--    surface fluxes if required. Note, each routine is called for different surface types. 
 !--    First call for default-type horizontal surfaces, for natural- and 
 !--    urban-type surfaces. Note, at this place only upward-facing horizontal
-!--    surfaces are treted. 
-       IF ( most_method == 'circular' )  THEN
-!
-!--       Default-type upward-facing horizontal surfaces
-          IF ( surf_def_h(0)%ns >= 1  )  THEN
-             surf => surf_def_h(0)
-             CALL calc_scaling_parameters
-             CALL calc_uvw_abs
-             IF ( .NOT. neutral )  CALL calc_ol
-             CALL calc_us
-             CALL calc_surface_fluxes
-             
-             IF ( do_output_at_2m )  THEN
-                CALL calc_pt_near_surface ( '2m' )
-             ENDIF
-          ENDIF
-!
-!--       Natural-type horizontal surfaces
-          IF ( surf_lsm_h%ns >= 1  )  THEN
-             surf => surf_lsm_h
-             CALL calc_scaling_parameters
-             CALL calc_uvw_abs
-             IF ( .NOT. neutral )  CALL calc_ol
-             CALL calc_us
-             CALL calc_surface_fluxes
-             IF ( do_output_at_2m )  THEN
-                CALL calc_pt_near_surface ( '2m' )
-             ENDIF
-          ENDIF
-!
-!--       Urban-type horizontal surfaces
-          IF ( surf_usm_h%ns >= 1  )  THEN
-             surf => surf_usm_h
-             CALL calc_scaling_parameters
-             CALL calc_uvw_abs
-             IF ( .NOT. neutral )  CALL calc_ol
-             CALL calc_us
-             CALL calc_surface_fluxes
-             IF ( do_output_at_2m )  THEN
-                CALL calc_pt_near_surface ( '2m' )
-             ENDIF
-             IF ( indoor_model )  THEN
-                CALL calc_pt_near_surface ( '10cm' )
-             ENDIF
-          ENDIF
-!
-!--    Use either Newton iteration or a lookup table for the bulk Richardson
-!--    number to calculate the Obukhov length 
-       ELSEIF ( most_method == 'newton'  .OR.  most_method == 'lookup' )  THEN
-!
-!--       Default-type upward-facing horizontal surfaces
-          IF ( surf_def_h(0)%ns >= 1  )  THEN
-             surf => surf_def_h(0)
-             CALL calc_uvw_abs
-             IF ( .NOT. neutral )  CALL calc_ol
-             CALL calc_us
-             CALL calc_scaling_parameters
-             CALL calc_surface_fluxes
-             IF ( do_output_at_2m )  THEN
-                CALL calc_pt_near_surface ( '2m' )
-             ENDIF
-          ENDIF
-!
-!--       Natural-type horizontal surfaces
-          IF ( surf_lsm_h%ns >= 1  )  THEN
-             surf => surf_lsm_h
-             CALL calc_uvw_abs
-             IF ( .NOT. neutral )  CALL calc_ol
-             CALL calc_us
-             CALL calc_scaling_parameters
-             CALL calc_surface_fluxes
-             IF ( do_output_at_2m )  THEN
-                CALL calc_pt_near_surface ( '2m' )
-             ENDIF
-          ENDIF
-!
-!--       Urban-type horizontal surfaces
-          IF ( surf_usm_h%ns >= 1  )  THEN
-             surf => surf_usm_h
-             CALL calc_uvw_abs
-             IF ( .NOT. neutral )  CALL calc_ol
-             CALL calc_us
-             CALL calc_scaling_parameters
-             CALL calc_surface_fluxes
-             IF ( do_output_at_2m )  THEN
-                CALL calc_pt_near_surface ( '2m' )
-             ENDIF
-             IF ( indoor_model )  THEN
-                CALL calc_pt_near_surface ( '10cm' )
-             ENDIF
-          ENDIF
-
+!--    surfaces are treated. 
+
+!
+!--    Default-type upward-facing horizontal surfaces
+       IF ( surf_def_h(0)%ns >= 1  )  THEN
+          surf => surf_def_h(0)
+          CALL calc_uvw_abs
+          IF ( .NOT. neutral )  CALL calc_ol
+          CALL calc_us
+          CALL calc_scaling_parameters
+          CALL calc_surface_fluxes
+          IF ( do_output_at_2m )  THEN
+             CALL calc_pt_near_surface ( '2m' )
+          ENDIF
        ENDIF
+!
+!--    Natural-type horizontal surfaces
+       IF ( surf_lsm_h%ns >= 1  )  THEN
+          surf => surf_lsm_h
+          CALL calc_uvw_abs
+          IF ( .NOT. neutral )  CALL calc_ol
+          CALL calc_us
+          CALL calc_scaling_parameters
+          CALL calc_surface_fluxes
+          IF ( do_output_at_2m )  THEN
+             CALL calc_pt_near_surface ( '2m' )
+          ENDIF
+       ENDIF
+!
+!--    Urban-type horizontal surfaces
+       IF ( surf_usm_h%ns >= 1  )  THEN
+          surf => surf_usm_h
+          CALL calc_uvw_abs
+          IF ( .NOT. neutral )  CALL calc_ol
+          CALL calc_us
+          CALL calc_scaling_parameters
+          CALL calc_surface_fluxes
+          IF ( do_output_at_2m )  THEN
+             CALL calc_pt_near_surface ( '2m' )
+          ENDIF
+          IF ( indoor_model )  THEN
+             CALL calc_pt_near_surface ( '10cm' )
+          ENDIF
+       ENDIF
+
 !
 !--    Treat downward-facing horizontal surfaces. Note, so far, these are 
@@ -549 +492 @@
 !--    flux in land-surface model in case of aero_resist_kray is not true.
        IF ( .NOT. aero_resist_kray )  THEN
-          IF ( most_method == 'circular' )  THEN
-             DO  l = 0, 1
-                IF ( surf_lsm_v(l)%ns >= 1  )  THEN
-                   surf => surf_lsm_v(l)
-!
-!--                Compute scaling parameters
-                   CALL calc_scaling_parameters
-!
-!--                Compute surface-parallel velocity
-                   CALL calc_uvw_abs_v_ugrid
-!
-!--                Compute Obukhov length
-                   IF ( .NOT. neutral )  CALL calc_ol
-!
-!--                Compute respective friction velocity on staggered grid
-                   CALL calc_us
-!
-!--                Compute respective surface fluxes for momentum and TKE
-                   CALL calc_surface_fluxes
-                ENDIF
-             ENDDO
-          ELSE
-             DO  l = 0, 1
-                IF ( surf_lsm_v(l)%ns >= 1  )  THEN
-                   surf => surf_lsm_v(l)
-!
-!--                Compute surface-parallel velocity
-                   CALL calc_uvw_abs_v_ugrid
-!
-!--                Compute Obukhov length
-                   IF ( .NOT. neutral )  CALL calc_ol
-!
-!--                Compute respective friction velocity on staggered grid
-                   CALL calc_us
-!
-!--                Compute scaling parameters
-                   CALL calc_scaling_parameters
-!
-!--                Compute respective surface fluxes for momentum and TKE
-                   CALL calc_surface_fluxes
-                ENDIF
-             ENDDO
-          ENDIF
+          DO  l = 0, 1
+             IF ( surf_lsm_v(l)%ns >= 1  )  THEN
+                surf => surf_lsm_v(l)
+!
+!--             Compute surface-parallel velocity
+                CALL calc_uvw_abs_v_ugrid
+!
+!--             Compute Obukhov length
+                IF ( .NOT. neutral )  CALL calc_ol
+!
+!--             Compute respective friction velocity on staggered grid
+                CALL calc_us
+!
+!--             Compute scaling parameters
+                CALL calc_scaling_parameters
+!
+!--             Compute respective surface fluxes for momentum and TKE
+                CALL calc_surface_fluxes
+             ENDIF
+          ENDDO
 !
 !--    No ts is required, so scaling parameters and Obukhov length do not need
@@ -652 +572 @@
 !--    flux in land-surface model in case of aero_resist_kray is not true.
        IF ( .NOT. aero_resist_kray )  THEN
-          IF ( most_method == 'circular' )  THEN
-             DO  l = 2, 3
-                IF ( surf_lsm_v(l)%ns >= 1  )  THEN
-                   surf => surf_lsm_v(l)
-!
-!--                Compute scaling parameters
-                   CALL calc_scaling_parameters
-!
-!--                Compute surface-parallel velocity
-                   CALL calc_uvw_abs_v_vgrid
-!
-!--                Compute Obukhov length
-                   IF ( .NOT. neutral )  CALL calc_ol
-!
-!--                Compute respective friction velocity on staggered grid
-                   CALL calc_us
-!
-!--                Compute respective surface fluxes for momentum and TKE
-                   CALL calc_surface_fluxes
-                ENDIF
-             ENDDO
-          ELSE
-             DO  l = 2, 3
-                IF ( surf_lsm_v(l)%ns >= 1  )  THEN
-                   surf => surf_lsm_v(l)
-!
-!--                Compute surface-parallel velocity
-                   CALL calc_uvw_abs_v_vgrid
-!
-!--                Compute Obukhov length
-                   IF ( .NOT. neutral )  CALL calc_ol
-!
-!--                Compute respective friction velocity on staggered grid
-                   CALL calc_us
-!
-!--                Compute scaling parameters
-                   CALL calc_scaling_parameters
-!
-!--                Compute respective surface fluxes for momentum and TKE
-                   CALL calc_surface_fluxes
-                ENDIF
-             ENDDO
-          ENDIF
+          DO  l = 2, 3
+             IF ( surf_lsm_v(l)%ns >= 1  )  THEN
+                surf => surf_lsm_v(l)
+!
+!--             Compute surface-parallel velocity
+                CALL calc_uvw_abs_v_vgrid
+!
+!--             Compute Obukhov length
+                IF ( .NOT. neutral )  CALL calc_ol
+!
+!--             Compute respective friction velocity on staggered grid
+                CALL calc_us
+!
+!--             Compute scaling parameters
+                CALL calc_scaling_parameters
+!
+!--             Compute respective surface fluxes for momentum and TKE
+                CALL calc_surface_fluxes
+             ENDIF
+          ENDDO
        ELSE
           DO  l = 2, 3
@@ -841 +738 @@
           IF ( surf_lsm_h%ns    >= 1 )  surf_lsm_h%ol    = 1.0E10_wp
           IF ( surf_usm_h%ns    >= 1 )  surf_usm_h%ol    = 1.0E10_wp
-       ENDIF
-
-       IF ( most_method == 'lookup' )  THEN
-
-!
-!--       Check for roughness heterogeneity. In that case terminate run and
-!--       inform user. Check for both, natural and non-natural walls.
-          IF ( surf_def_h(0)%ns >= 1 )  THEN
-             IF ( MINVAL( surf_def_h(0)%z0h ) /= MAXVAL( surf_def_h(0)%z0h )  .OR. &
-                  MINVAL( surf_def_h(0)%z0  ) /= MAXVAL( surf_def_h(0)%z0  ) )  THEN
-                terminate_run_l = .TRUE.
-             ENDIF
-          ENDIF
-          IF ( surf_lsm_h%ns >= 1 )  THEN
-             IF ( MINVAL( surf_lsm_h%z0h ) /= MAXVAL( surf_lsm_h%z0h )  .OR.       &
-                  MINVAL( surf_lsm_h%z0  ) /= MAXVAL( surf_lsm_h%z0  ) )  THEN
-                terminate_run_l = .TRUE.
-             ENDIF
-          ENDIF
-          IF ( surf_usm_h%ns >= 1 )  THEN
-             IF ( MINVAL( surf_usm_h%z0h ) /= MAXVAL( surf_usm_h%z0h )  .OR.       &
-                  MINVAL( surf_usm_h%z0  ) /= MAXVAL( surf_usm_h%z0  ) )  THEN
-                terminate_run_l = .TRUE.
-             ENDIF
-          ENDIF
-!
-!--       Check roughness homogeneity between differt surface types. 
-          IF ( surf_lsm_h%ns >= 1  .AND.  surf_def_h(0)%ns >= 1 )  THEN
-             IF ( MINVAL( surf_lsm_h%z0h ) /= MAXVAL( surf_def_h(0)%z0h )  .OR.    &
-                  MINVAL( surf_lsm_h%z0  ) /= MAXVAL( surf_def_h(0)%z0  ) )  THEN
-                terminate_run_l = .TRUE.
-             ENDIF
-          ENDIF
-          IF ( surf_usm_h%ns >= 1  .AND.  surf_def_h(0)%ns >= 1 )  THEN
-             IF ( MINVAL( surf_usm_h%z0h ) /= MAXVAL( surf_def_h(0)%z0h )  .OR.    &
-                  MINVAL( surf_usm_h%z0  ) /= MAXVAL( surf_def_h(0)%z0  ) )  THEN
-                terminate_run_l = .TRUE.
-             ENDIF
-          ENDIF
-          IF ( surf_usm_h%ns >= 1  .AND.  surf_lsm_h%ns >= 1 )  THEN
-             IF ( MINVAL( surf_usm_h%z0h ) /= MAXVAL( surf_lsm_h%z0h )  .OR.       &
-                  MINVAL( surf_usm_h%z0  ) /= MAXVAL( surf_lsm_h%z0  ) )  THEN
-                terminate_run_l = .TRUE.
-             ENDIF
-          ENDIF
-
-
-#if defined( __parallel )
-!
-!--       Make a logical OR for all processes. Force termiation of model if result
-!--       is TRUE
-          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
-          CALL MPI_ALLREDUCE( terminate_run_l, terminate_run, 1, MPI_LOGICAL,  &
-                              MPI_LOR, comm2d, ierr )
-#else
-          terminate_run = terminate_run_l
-#endif
-
-          IF ( terminate_run )  THEN
-             message_string = 'most_method = "lookup" cannot be used in ' //   &
-                              'combination with a prescribed roughness '  //   &
-                              'heterogeneity'
-             CALL message( 'surface_layer_fluxes', 'PA0116', 1, 2, 0, 6, 0 )
-          ENDIF
-
-          ALLOCATE(  zeta_tmp(0:num_steps-1) )
-
-!
-!--       Use the lowest possible value for z_mo
-          k    = nzb
-          z_mo = zu(k+1) - zw(k)
-
-!
-!--       Calculate z/L range from zeta_stretch to zeta_max using 90% of the
-!--       available steps (num_steps). The calculation is done with negative
-!--       values of zeta in order to simplify the stretching in the free
-!--       convection limit for the remaining 10% of steps.
-          zeta_tmp(0) = - zeta_max
-          num_steps_n = ( num_steps * 9 / 10 ) - 1
-          zeta_step   = (zeta_max - zeta_stretch) / REAL(num_steps_n)
-
-          DO li = 1, num_steps_n
-             zeta_tmp(li) = zeta_tmp(li-1) + zeta_step
-          ENDDO
-
-!
-!--       Calculate stretching factor for the free convection range
-          DO  WHILE ( ABS( (regr-regr_old) / regr_old ) > 1.0E-10_wp )
-             regr_old = regr
-             regr = ( 1.0_wp - ( -zeta_min / zeta_step ) * ( 1.0_wp - regr )   &
-                    )**( 10.0_wp / REAL(num_steps) )
-          ENDDO
-
-!
-!--       Calculate z/L range from zeta_min to zeta_stretch
-          DO li = num_steps_n+1, num_steps-1
-             zeta_tmp(li) = zeta_tmp(li-1) + zeta_step
-             zeta_step = zeta_step * regr
-          ENDDO
-
-!
-!--       Save roughness lengths to temporary variables
-          IF ( surf_def_h(0)%ns >= 1  )  THEN
-             z0h_min = surf_def_h(0)%z0h(1)
-             z0_min  = surf_def_h(0)%z0(1)
-          ELSEIF ( surf_lsm_h%ns >= 1 )  THEN
-             z0h_min = surf_lsm_h%z0h(1)
-             z0_min  = surf_lsm_h%z0(1)
-          ELSEIF ( surf_usm_h%ns >= 1 )  THEN
-             z0h_min = surf_usm_h%z0h(1)
-             z0_min  = surf_usm_h%z0(1)
-          ENDIF         
-!
-!--       Calculate lookup table for the Richardson number versus Obukhov length
-!--       The Richardson number (rib) is defined depending on the choice of 
-!--       boundary conditions for temperature
-          IF ( ibc_pt_b == 1 )  THEN
-             DO li = 0, num_steps-1
-                ol_tab(li)  = - z_mo / zeta_tmp(num_steps-1-li)
-                rib_tab(li) = z_mo / ol_tab(li)  / ( LOG( z_mo / z0_min )       &
-                                                - psi_m( z_mo / ol_tab(li) )    &
-                                                + psi_m( z0_min / ol_tab(li) )  &
-                                                  )**3
-             ENDDO  
-          ELSE
-             DO li = 0, num_steps-1
-                ol_tab(li)  = - z_mo / zeta_tmp(num_steps-1-li)
-                rib_tab(li) = z_mo / ol_tab(li)  * ( LOG( z_mo / z0h_min )     &
-                                              - psi_h( z_mo / ol_tab(li) )     &
-                                              + psi_h( z0h_min / ol_tab(li) )  &
-                                            )                                  &
-                                          / ( LOG( z_mo / z0_min )             &
-                                              - psi_m( z_mo / ol_tab(li) )     &
-                                              + psi_m( z0_min / ol_tab(li) )   &
-                                            )**2
-             ENDDO
-          ENDIF
-
-!
-!--       Determine minimum values of rib in the lookup table. Set upper limit 
-!--       to critical Richardson number (0.25)
-          rib_min  = MINVAL(rib_tab)
-          rib_max  = 0.25 !MAXVAL(rib_tab)
-
-          DEALLOCATE( zeta_tmp )
        ENDIF
 
@@ -1256 +1008 @@
                        ol_u      !< Upper bound of L for Newton iteration
 
-       IF ( TRIM( most_method ) /= 'circular' )  THEN
-!
-!--       Evaluate bulk Richardson number (calculation depends on
-!--       definition based on setting of boundary conditions
-          IF ( ibc_pt_b /= 1 )  THEN
-             IF ( humidity )  THEN
-                !$OMP PARALLEL DO PRIVATE( z_mo )
-                DO  m = 1, surf%ns
-
-                   z_mo = surf%z_mo(m)
-
-                   surf%rib(m) = g * z_mo                                      &
-                                 * ( surf%vpt1(m) - surf%vpt_surface(m) )      &
-                                 / ( surf%uvw_abs(m)**2 * surf%vpt1(m)         &
-                                 + 1.0E-20_wp )
-                ENDDO
-             ELSE
-                !$OMP PARALLEL DO PRIVATE( z_mo )
-                DO  m = 1, surf%ns
-
-                   z_mo = surf%z_mo(m)
-
-                   surf%rib(m) = g * z_mo                                      &
-                                 * ( surf%pt1(m) - surf%pt_surface(m) )        &
-                                 / ( surf%uvw_abs(m)**2 * surf%pt1(m)          &
-                                 + 1.0E-20_wp )
-                ENDDO
-             ENDIF
+!
+!--    Evaluate bulk Richardson number (calculation depends on
+!--    definition based on setting of boundary conditions
+       IF ( ibc_pt_b /= 1 )  THEN
+          IF ( humidity )  THEN
+             !$OMP PARALLEL DO PRIVATE( z_mo )
+             DO  m = 1, surf%ns
+
+                z_mo = surf%z_mo(m)
+
+                surf%rib(m) = g * z_mo                                         &
+                              * ( surf%vpt1(m) - surf%vpt_surface(m) )         &
+                              / ( surf%uvw_abs(m)**2 * surf%vpt1(m)            &
+                              + 1.0E-20_wp )
+             ENDDO
           ELSE
-             IF ( humidity )  THEN
-                !$OMP PARALLEL DO PRIVATE( k, z_mo )
-                DO  m = 1, surf%ns
-
-                   k   = surf%k(m)
-
-                   z_mo = surf%z_mo(m)
-
-                   surf%rib(m) = - g * z_mo * ( ( 1.0_wp + 0.61_wp             &
+             !$OMP PARALLEL DO PRIVATE( z_mo )
+             DO  m = 1, surf%ns
+
+                z_mo = surf%z_mo(m)
+
+                surf%rib(m) = g * z_mo                                         &
+                              * ( surf%pt1(m) - surf%pt_surface(m) )           &
+                              / ( surf%uvw_abs(m)**2 * surf%pt1(m) + 1.0E-20_wp )
+             ENDDO
+          ENDIF
+       ELSE
+          IF ( humidity )  THEN
+             !$OMP PARALLEL DO PRIVATE( k, z_mo )
+             DO  m = 1, surf%ns
+
+                k   = surf%k(m)
+
+                z_mo = surf%z_mo(m)
+
+                surf%rib(m) = - g * z_mo * ( ( 1.0_wp + 0.61_wp                &
                            * surf%qv1(m) ) * surf%shf(m) + 0.61_wp             &
                            * surf%pt1(m) * surf%qsws(m) ) *                    &
@@ -1299 +1049 @@
                          ( surf%uvw_abs(m)**3 * surf%vpt1(m) * kappa**2        &
                            + 1.0E-20_wp )
-                ENDDO
+             ENDDO
+          ELSE
+             !$OMP PARALLEL DO PRIVATE( k, z_mo )
+             !$ACC PARALLEL LOOP PRIVATE(k, z_mo) &
+             !$ACC PRESENT(surf, drho_air_zw)
+             DO  m = 1, surf%ns
+
+                k   = surf%k(m)
+
+                z_mo = surf%z_mo(m)
+
+                surf%rib(m) = - g * z_mo * surf%shf(m) *                    &
+                                     drho_air_zw(k-1)            /          &
+                        ( surf%uvw_abs(m)**3 * surf%pt1(m) * kappa**2       &
+                        + 1.0E-20_wp )
+             ENDDO
+          ENDIF
+       ENDIF
+
+!
+!--    Calculate the Obukhov length using Newton iteration
+       !$ACC PARALLEL LOOP PRIVATE(i, j, z_mo) &
+       !$ACC PRIVATE(ol_old, ol_m, ol_l, ol_u, f, f_d_ol) &
+       !$ACC PRESENT(surf)
+       DO  m = 1, surf%ns
+
+          i   = surf%i(m)            
+          j   = surf%j(m)
+
+          z_mo = surf%z_mo(m)
+
+!
+!--       Store current value in case the Newton iteration fails
+          ol_old = surf%ol(m)
+
+!
+!--       Ensure that the bulk Richardson number and the Obukhov 
+!--       length have the same sign
+          IF ( surf%rib(m) * surf%ol(m) < 0.0_wp  .OR.                      &
+               ABS( surf%ol(m) ) == ol_max )  THEN
+             IF ( surf%rib(m) > 1.0_wp ) surf%ol(m) =  0.01_wp
+             IF ( surf%rib(m) < 0.0_wp ) surf%ol(m) = -0.01_wp
+          ENDIF
+!
+!--       Iteration to find Obukhov length
+          iter = 0
+          DO
+             iter = iter + 1
+!
+!--          In case of divergence, use the value of the previous time step
+             IF ( iter > 1000 )  THEN
+                surf%ol(m) = ol_old
+                EXIT
+             ENDIF
+
+             ol_m = surf%ol(m)
+             ol_l = ol_m - 0.001_wp * ol_m
+             ol_u = ol_m + 0.001_wp * ol_m
+
+
+             IF ( ibc_pt_b /= 1 )  THEN
+!
+!--             Calculate f = Ri - [...]/[...]^2 = 0
+                f = surf%rib(m) - ( z_mo / ol_m ) * (                          &
+                                              LOG( z_mo / surf%z0h(m) )        &
+                                              - psi_h( z_mo / ol_m )           &
+                                              + psi_h( surf%z0h(m) /           &
+                                                       ol_m )                  &  
+                                                     )                         &
+                                           / ( LOG( z_mo / surf%z0(m) )        &
+                                              - psi_m( z_mo / ol_m )           &
+                                              + psi_m( surf%z0(m) /  ol_m )    &
+                                                 )**2
+
+!
+!--              Calculate df/dL
+                 f_d_ol = ( - ( z_mo / ol_u ) * ( LOG( z_mo /                  &
+                                                          surf%z0h(m) )        &
+                                         - psi_h( z_mo / ol_u )                &
+                                         + psi_h( surf%z0h(m) / ol_u )         &
+                                           )                                   &
+                                         / ( LOG( z_mo / surf%z0(m) )          &
+                                         - psi_m( z_mo / ol_u )                &
+                                         + psi_m( surf%z0(m) / ol_u )          &
+                                           )**2                                &
+                        + ( z_mo / ol_l ) * ( LOG( z_mo / surf%z0h(m) )        &
+                                         - psi_h( z_mo / ol_l )                &
+                                         + psi_h( surf%z0h(m) / ol_l )         &
+                                           )                                   &
+                                         / ( LOG( z_mo / surf%z0(m) )          &
+                                         - psi_m( z_mo / ol_l )                &
+                                         + psi_m( surf%z0(m) / ol_l )          &
+                                           )**2                                &
+                          ) / ( ol_u - ol_l )
              ELSE
-                !$OMP PARALLEL DO PRIVATE( k, z_mo )
-                !$ACC PARALLEL LOOP PRIVATE(k, z_mo) &
-                !$ACC PRESENT(surf, drho_air_zw)
-                DO  m = 1, surf%ns
-
-                   k   = surf%k(m)
-
-                   z_mo = surf%z_mo(m)
-
-                   surf%rib(m) = - g * z_mo * surf%shf(m) *                    &
-                                        drho_air_zw(k-1)            /          &
-                        ( surf%uvw_abs(m)**3 * surf%pt1(m) * kappa**2          &
-                           + 1.0E-20_wp )
-                ENDDO
+!
+!--             Calculate f = Ri - 1 /[...]^3 = 0
+                f = surf%rib(m) - ( z_mo / ol_m ) /                            &
+                                             ( LOG( z_mo / surf%z0(m) )        &
+                                         - psi_m( z_mo / ol_m )                &
+                                         + psi_m( surf%z0(m) / ol_m )          &
+                                             )**3
+
+!
+!--             Calculate df/dL
+                f_d_ol = ( - ( z_mo / ol_u ) / ( LOG( z_mo / surf%z0(m) )      &
+                                         - psi_m( z_mo / ol_u )                &
+                                         + psi_m( surf%z0(m) / ol_u )          &
+                                                  )**3                         &
+                           + ( z_mo / ol_l ) / ( LOG( z_mo / surf%z0(m) )      &
+                                         - psi_m( z_mo / ol_l )                &
+                                         + psi_m( surf%z0(m) / ol_l )          &
+                                            )**3                               &
+                          ) / ( ol_u - ol_l )
              ENDIF
-          ENDIF
-
-       ENDIF
-
-
-!
-!--    Calculate the Obukhov length either using a Newton iteration
-!--    method, via a lookup table, or using the old circular way
-       IF ( TRIM( most_method ) == 'newton' )  THEN
-
-          !$ACC PARALLEL LOOP PRIVATE(i, j, z_mo) &
-          !$ACC PRIVATE(ol_old, ol_m, ol_l, ol_u, f, f_d_ol) &
-          !$ACC PRESENT(surf)
-          DO  m = 1, surf%ns
-
-             i   = surf%i(m)            
-             j   = surf%j(m)
-
-             z_mo = surf%z_mo(m)
-
-!
-!--          Store current value in case the Newton iteration fails
-             ol_old = surf%ol(m)
+!
+!--          Calculate new L
+             surf%ol(m) = ol_m - f / f_d_ol
 
 !
 !--          Ensure that the bulk Richardson number and the Obukhov 
-!--          length have the same sign
-             IF ( surf%rib(m) * surf%ol(m) < 0.0_wp  .OR.                      &
-                  ABS( surf%ol(m) ) == ol_max )  THEN
-                IF ( surf%rib(m) > 1.0_wp ) surf%ol(m) =  0.01_wp
-                IF ( surf%rib(m) < 0.0_wp ) surf%ol(m) = -0.01_wp
+!--          length have the same sign and ensure convergence.
+             IF ( surf%ol(m) * ol_m < 0.0_wp )  surf%ol(m) = ol_m * 0.5_wp
+!
+!--          If unrealistic value occurs, set L to the maximum
+!--          value that is allowed
+             IF ( ABS( surf%ol(m) ) > ol_max )  THEN
+                surf%ol(m) = ol_max
+                EXIT
              ENDIF
 !
-!--          Iteration to find Obukhov length
-             iter = 0
-             DO
-                iter = iter + 1
-!
-!--             In case of divergence, use the value of the previous time step
-                IF ( iter > 1000 )  THEN
-                   surf%ol(m) = ol_old
-                   EXIT
-                ENDIF
-
-                ol_m = surf%ol(m)
-                ol_l = ol_m - 0.001_wp * ol_m
-                ol_u = ol_m + 0.001_wp * ol_m
-
-
-                IF ( ibc_pt_b /= 1 )  THEN
-!
-!--                Calculate f = Ri - [...]/[...]^2 = 0
-                   f = surf%rib(m) - ( z_mo / ol_m ) * (                       &
-                                                 LOG( z_mo / surf%z0h(m) )     &
-                                                 - psi_h( z_mo / ol_m )        &
-                                                 + psi_h( surf%z0h(m) /        &
-                                                          ol_m )               &
-                                                               )               &
-                                              / ( LOG( z_mo / surf%z0(m) )     &
-                                                 - psi_m( z_mo / ol_m )        &
-                                                 + psi_m( surf%z0(m) /         &
-                                                          ol_m )               &
-                                                 )**2
-
-!
-!--                 Calculate df/dL
-                    f_d_ol = ( - ( z_mo / ol_u ) * ( LOG( z_mo /               &
-                                                             surf%z0h(m) )     &
-                                            - psi_h( z_mo / ol_u )             &
-                                            + psi_h( surf%z0h(m) / ol_u )      &
-                                              )                                &
-                                            / ( LOG( z_mo / surf%z0(m) )       &
-                                            - psi_m( z_mo / ol_u )             &
-                                            + psi_m( surf%z0(m) / ol_u )       &
-                                              )**2                             &
-                           + ( z_mo / ol_l ) * ( LOG( z_mo / surf%z0h(m) )     &
-                                            - psi_h( z_mo / ol_l )             &
-                                            + psi_h( surf%z0h(m) / ol_l )      &
-                                              )                                &
-                                            / ( LOG( z_mo / surf%z0(m) )       &
-                                            - psi_m( z_mo / ol_l )             &
-                                            + psi_m( surf%z0(m) / ol_l )       &
-                                              )**2                             &
-                             ) / ( ol_u - ol_l )
-                ELSE
-!
-!--                Calculate f = Ri - 1 /[...]^3 = 0
-                   f = surf%rib(m) - ( z_mo / ol_m ) /                         &
-                                                ( LOG( z_mo / surf%z0(m) )     &
-                                            - psi_m( z_mo / ol_m )             &
-                                            + psi_m( surf%z0(m) / ol_m )       &
-                                                )**3
-
-!
-!--                Calculate df/dL
-                   f_d_ol = ( - ( z_mo / ol_u ) / ( LOG( z_mo /                &
-                                                            surf%z0(m) )       &
-                                            - psi_m( z_mo / ol_u )             &
-                                            + psi_m( surf%z0(m) / ol_u )       &
-                                                     )**3                      &
-                           + ( z_mo / ol_l ) / ( LOG( z_mo / surf%z0(m) )      &
-                                            - psi_m( z_mo / ol_l )             &
-                                            + psi_m( surf%z0(m) / ol_l )       &
-                                               )**3                            &
-                                  ) / ( ol_u - ol_l )
-                ENDIF
-!
-!--             Calculate new L
-                surf%ol(m) = ol_m - f / f_d_ol
-
-!
-!--             Ensure that the bulk Richardson number and the Obukhov 
-!--             length have the same sign and ensure convergence.
-                IF ( surf%ol(m) * ol_m < 0.0_wp )  surf%ol(m) = ol_m * 0.5_wp
-!
-!--             If unrealistic value occurs, set L to the maximum
-!--             value that is allowed
-                IF ( ABS( surf%ol(m) ) > ol_max )  THEN
-                   surf%ol(m) = ol_max
-                   EXIT
-                ENDIF
-!
-!--             Check for convergence
-                IF ( ABS( ( surf%ol(m) - ol_m ) /                              &
-                     surf%ol(m) ) < 1.0E-4_wp )  THEN
-                   EXIT
-                ELSE
-                   CYCLE
-                ENDIF
-
-             ENDDO
+!--          Check for convergence
+             IF ( ABS( ( surf%ol(m) - ol_m ) /  surf%ol(m) ) < 1.0E-4_wp )  THEN
+                EXIT
+             ELSE
+                CYCLE
+             ENDIF
+
           ENDDO
-
-       ELSEIF ( TRIM( most_method ) == 'lookup' )  THEN
-
-          !$OMP PARALLEL DO PRIVATE( i, j, z_mo, li ) FIRSTPRIVATE( li_bnd ) LASTPRIVATE( li_bnd )
-          DO  m = 1, surf%ns
-
-             i   = surf%i(m)            
-             j   = surf%j(m)
-!
-!--          If the bulk Richardson number is outside the range of the lookup
-!--          table, set it to the exceeding threshold value
-             IF ( surf%rib(m) < rib_min )  surf%rib(m) = rib_min
-             IF ( surf%rib(m) > rib_max )  surf%rib(m) = rib_max
-
-!
-!--          Find the correct index bounds for linear interpolation. As the
-!--          Richardson number will not differ very much from time step to
-!--          time step , use the index from the last step and search in the 
-!--          correct direction
-             li = li_bnd
-             IF ( rib_tab(li) - surf%rib(m) > 0.0_wp )  THEN
-                DO  WHILE ( rib_tab(li-1) - surf%rib(m) > 0.0_wp  .AND.  li > 0 )
-                   li = li-1
-                ENDDO
-             ELSE
-                DO  WHILE ( rib_tab(li) - surf%rib(m) < 0.0_wp                 &
-                           .AND.  li < num_steps-1 )
-                   li = li+1
-                ENDDO
-             ENDIF
-             li_bnd = li
-
-!
-!--          Linear interpolation to find the correct value of z/L
-             surf%ol(m) = ( ol_tab(li-1) + ( ol_tab(li) - ol_tab(li-1) )       &
-                         / (  rib_tab(li) - rib_tab(li-1) )                    &
-                         * ( surf%rib(m) - rib_tab(li-1) ) )
-
-          ENDDO
-
-       ELSEIF ( TRIM( most_method ) == 'circular' )  THEN
-
-          IF ( .NOT. humidity )  THEN
-             !$OMP PARALLEL DO PRIVATE( z_mo )
-             DO  m = 1, surf%ns
-
-                z_mo = surf%z_mo(m)
-
-                surf%ol(m) =  ( surf%pt1(m) *  surf%us(m)**2 ) /               &
-                                       ( kappa * g *                           &
-                                         surf%ts(m) + 1E-30_wp )
-!
-!--             Limit the value range of the Obukhov length.
-!--             This is necessary for very small velocities (u,v --> 1), because
-!--             the absolute value of ol can then become very small, which in
-!--             consequence would result in very large shear stresses and very
-!--             small momentum fluxes (both are generally unrealistic).
-                IF ( ( z_mo / ( surf%ol(m) + 1E-30_wp ) ) < zeta_min )         &
-                   surf%ol(m) = z_mo / zeta_min
-                IF ( ( z_mo / ( surf%ol(m) + 1E-30_wp ) ) > zeta_max )         &
-                   surf%ol(m) = z_mo / zeta_max
-
-             ENDDO
-          ELSEIF ( bulk_cloud_model  .OR.  cloud_droplets )  THEN
-             !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
-             DO  m = 1, surf%ns
-
-                i   = surf%i(m)            
-                j   = surf%j(m)
-                k   = surf%k(m)
-
-                z_mo = surf%z_mo(m)
-
-
-                surf%ol(m) =  ( surf%vpt1(m) * surf%us(m)**2 ) /               &
-                    ( kappa * g * ( surf%ts(m) +                               &
-                         0.61_wp * surf%pt1(m) * surf%us(m)                    &
-                       + 0.61_wp * surf%qv1(m) * surf%ts(m) -                  &
-                                   surf%ts(m)  * ql(k,j,i) ) + 1E-30_wp )
-!
-!--             Limit the value range of the Obukhov length.
-!--             This is necessary for very small velocities (u,v --> 1), because
-!--             the absolute value of ol can then become very small, which in
-!--             consequence would result in very large shear stresses and very
-!--             small momentum fluxes (both are generally unrealistic).
-                IF ( ( z_mo / ( surf%ol(m) + 1E-30_wp ) ) < zeta_min )         &
-                   surf%ol(m) = z_mo / zeta_min
-                IF ( ( z_mo / ( surf%ol(m) + 1E-30_wp ) ) > zeta_max )         &
-                   surf%ol(m) = z_mo / zeta_max
-
-             ENDDO
-          ELSE
-
-             !$OMP PARALLEL DO PRIVATE( z_mo )
-             DO  m = 1, surf%ns
-
-                z_mo = surf%z_mo(m)
-
-                surf%ol(m) =  ( surf%vpt1(m) *  surf%us(m)**2 ) /              &
-                    ( kappa * g * ( surf%ts(m) + 0.61_wp * surf%pt1(m) *       &
-                                    surf%qs(m) + 0.61_wp * surf%qv1(m) *       &
-                                    surf%ts(m) ) + 1E-30_wp )
-
-!
-!--             Limit the value range of the Obukhov length.
-!--             This is necessary for very small velocities (u,v --> 1), because
-!--             the absolute value of ol can then become very small, which in
-!--             consequence would result in very large shear stresses and very
-!--             small momentum fluxes (both are generally unrealistic).
-                IF ( ( z_mo / ( surf%ol(m) + 1E-30_wp ) ) < zeta_min )         &
-                   surf%ol(m) = z_mo / zeta_min
-                IF ( ( z_mo / ( surf%ol(m) + 1E-30_wp ) ) > zeta_max )         &
-                   surf%ol(m) = z_mo / zeta_max
-
-             ENDDO
-
-          ENDIF
-
-       ENDIF
+       ENDDO
 
     END SUBROUTINE calc_ol

palm/trunk/SOURCE/write_restart_data_mod.f90

@@ -163 +163 @@
 
 
-       binary_version_global = '4.7'
+       binary_version_global = '4.8'
 
        CALL wrd_write_string( 'binary_version_global' )
@@ -466 +466 @@
        WRITE ( 14 )  momentum_advec
 
-       CALL wrd_write_string( 'most_method' ) 
-       WRITE ( 14 )  most_method
-
        CALL wrd_write_string( 'netcdf_precision' ) 
        WRITE ( 14 )  netcdf_precision