Changeset 3668


Ignore:
Timestamp:
Jan 14, 2019 12:49:24 PM (3 years ago)
Author:
maronga
Message:

removed most_methods circular and lookup. added improved version of palm_csd

Location:
palm/trunk
Files:
12 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SCRIPTS/.csd.config.default

    r3567 r3668  
    2525# -----------------
    2626# $Id$
     27# Updated configuration
     28#
     29# 3567 2018-11-27 13:59:21Z maronga
    2730# Initial revisions
    28 #
    29 #
    30 #
    31 #
    3231#
    3332# Description:
     
    5251institution = Institute of Meteorology and Climatology, Leibniz University Hannover
    5352palm_version = 6.0
    54 version = 1
    5553rotation_angle = 0.0
    5654references = text
     
    5856
    5957[settings]
    60 filename_out = /ldata2/MOSAIK/showcase_python
    61 lai_season = spring
    62 lai_high_vegetation_default = 5.0
    63 lai_low_vegetation_default = 1.0
     58lai_roof_intensive = 1.5
     59lai_roof_extensive = 3.0
     60lai_high_vegetation_default = 6.0
     61lai_low_vegetation_default = 3.0
    6462lai_alpha = 5.0
    6563lai_beta = 3.0
    66 tree_height_default = 10.0
     64patch_height_default = 10.0
    6765bridge_width = 3.0
    6866debug_mode = False
    69 
     67season = "summer"
    7068
    7169[input_01]
     
    7977file_lon = Berlin_CoordinatesLatLon_x_15m_DLR.nc
    8078file_zt = Berlin_terrain_height_15m_DLR.nc
    81 file_buildings_2d = Berlin_building_height_15m_LOD2_DLR.nc
     79file_buildings_2d = Berlin_building_height_15m_DLR.nc
    8280file_building_id = Berlin_building_id_15m_DLR.nc
    8381file_building_type = Berlin_building_type_15m_DLR.nc
    8482file_bridges_2d = Berlin_bridges_height_15m_DLR.nc
    8583file_bridges_id = Berlin_bridges_id_15m_DLR.nc
     84file_lai =  Berlin_leaf_area_index_15m_DLR_WANG_summer.nc
    8685file_vegetation_type = Berlin_vegetation_type_15m_DLR.nc
    8786file_vegetation_height = Berlin_vegetation_patch_height_15m_DLR.nc
     
    9190file_street_type =  Berlin_street_type_15m_DLR.nc
    9291file_street_crossings =  Berlin_street_crossings_15m_DLR.nc
     92file_tree_height =  Berlin_trees_height_clean_15m_DLR.nc
     93file_tree_crown_diameter =  Berlin_trees_crown_clean_15m_DLR.nc
     94file_tree_trunk_diameter =  Berlin_trees_trunk_clean_15m.nc
     95file_tree_type =  Berlin_trees_type_15m_DLR.nc
     96file_patch_height =  Berlin_vegetation_patch_height_15m_DLR.nc
     97file_vegetation_on_roofs =  Berlin_vegetation_on_roofs_15m_DLR.nc
    9398
    9499[input_02]
     
    102107file_lon = Berlin_CoordinatesLatLon_x_10m_DLR.nc
    103108file_zt = Berlin_terrain_height_10m_DLR.nc
    104 file_buildings_2d = Berlin_building_height_10m_LOD2_DLR.nc
     109file_buildings_2d = Berlin_building_height_10m_DLR.nc
    105110file_building_id = Berlin_building_id_10m_DLR.nc
    106111file_building_type = Berlin_building_type_10m_DLR.nc
    107112file_bridges_2d = Berlin_bridges_height_10m_DLR.nc
    108113file_bridges_id = Berlin_bridges_id_10m_DLR.nc
     114file_lai =  Berlin_leaf_area_index_10m_DLR_WANG_summer.nc
    109115file_vegetation_type = Berlin_vegetation_type_10m_DLR.nc
    110116file_vegetation_height = Berlin_vegetation_patch_height_10m_DLR.nc
     
    114120file_street_type =  Berlin_street_type_10m_DLR.nc
    115121file_street_crossings =  Berlin_street_crossings_10m_DLR.nc
     122file_tree_height =  Berlin_trees_height_clean_10m_DLR.nc
     123file_tree_crown_diameter =  Berlin_trees_crown_10m_DLR.nc
     124file_tree_trunk_diameter =  Berlin_trees_trunk_clean_10m.nc
     125file_tree_type =  Berlin_trees_type_10m_DLR.nc
     126file_patch_height =  Berlin_vegetation_patch_height_10m_DLR.nc
     127file_vegetation_on_roofs =  Berlin_vegetation_on_roofs_10m_DLR.nc
    116128
    117129[input_03]
     130path = /ldata2/MOSAIK/Berlin_static_driver_data
     131pixel_size = 1.0
    118132file_x = Berlin_CoordinatesUTM_x_1m_DLR.nc
    119133file_y = Berlin_CoordinatesUTM_y_1m_DLR.nc
     
    122136file_lat = Berlin_CoordinatesLatLon_y_1m_DLR.nc
    123137file_lon = Berlin_CoordinatesLatLon_x_1m_DLR.nc
    124 path = /ldata2/MOSAIK/Berlin_static_driver_data
    125 pixel_size = 1.0
    126138file_zt = Berlin_terrain_height_1m_DLR.nc
    127 file_buildings_2d = Berlin_building_height_1m_LOD2_DLR.nc
     139file_buildings_2d = Berlin_building_height_1m_DLR.nc
    128140file_building_id = Berlin_building_id_1m_DLR.nc
    129141file_bridges_2d = Berlin_bridges_height_1m_DLR.nc
    130142file_bridges_id = Berlin_bridges_id_1m_DLR.nc
    131143file_building_type = Berlin_building_type_1m_DLR.nc
     144file_lai =  Berlin_leaf_area_index_1m_DLR_WANG_summer.nc
    132145file_vegetation_type = Berlin_vegetation_type_1m_DLR.nc
    133146file_vegetation_height = Berlin_vegetation_patch_height_1m_DLR.nc
     
    137150file_street_type =  Berlin_street_type_1m_DLR.nc
    138151file_street_crossings =  Berlin_street_crossings_1m_DLR.nc
     152file_tree_height =  Berlin_trees_height_clean_1m_DLR.nc
     153file_tree_crown_diameter =  Berlin_trees_crown_clean_1m_DLR.nc
     154file_tree_trunk_diameter =  Berlin_trees_trunk_clean_1m.nc
     155file_tree_type =  Berlin_trees_type_1m_DLR.nc
     156file_patch_height =  Berlin_vegetation_patch_height_1m_DLR.nc
     157file_vegetation_on_roofs =  Berlin_vegetation_on_roofs_1m_DLR.nc
     158
     159[input_04]
     160path = /ldata2/MOSAIK/Berlin_static_driver_data
     161pixel_size = 2.0
     162file_x = Berlin_CoordinatesUTM_x_2m_DLR.nc
     163file_y = Berlin_CoordinatesUTM_y_2m_DLR.nc
     164file_x_UTM = Berlin_CoordinatesUTM_x_2m_DLR.nc
     165file_y_UTM = Berlin_CoordinatesUTM_y_2m_DLR.nc
     166file_lat = Berlin_CoordinatesLatLon_y_2m_DLR.nc
     167file_lon = Berlin_CoordinatesLatLon_x_2m_DLR.nc
     168file_zt = Berlin_terrain_height_2m_DLR.nc
     169file_buildings_2d = Berlin_building_height_2m_DLR.nc
     170file_building_id = Berlin_building_id_2m_DLR.nc
     171file_bridges_2d = Berlin_bridges_height_2m_DLR.nc
     172file_bridges_id = Berlin_bridges_id_2m_DLR.nc
     173file_building_type = Berlin_building_type_2m_DLR.nc
     174file_lai =  Berlin_leaf_area_index_2m_DLR_WANG_summer.nc
     175file_vegetation_type = Berlin_vegetation_type_2m_DLR.nc
     176file_vegetation_height = Berlin_vegetation_patch_height_2m_DLR.nc
     177file_pavement_type = Berlin_pavement_type_2m_DLR.nc
     178file_water_type = Berlin_water_type_2m_DLR.nc
     179file_soil_type = Berlin_soil_type_2m_DLR.nc
     180file_street_type =  Berlin_street_type_2m_DLR.nc
     181file_street_crossings =  Berlin_street_crossings_2m_DLR.nc
     182file_tree_height =  Berlin_trees_height_clean_2m_DLR.nc
     183file_tree_crown_diameter = Berlin_trees_crown_clean_2m_DLR.nc
     184file_tree_trunk_diameter =  Berlin_trees_trunk_clean_2m.nc
     185file_tree_type =  Berlin_trees_type_2m_DLR.nc
     186file_patch_height =  Berlin_vegetation_patch_height_2m_DLR.nc
     187file_vegetation_on_roofs =  Berlin_vegetation_on_roofs_2m_DLR.nc
     188
    139189
    140190[output]
    141191path = /ldata2/MOSAIK/
    142 file_out = winter_iop1_test
     192file_out = showcase_berlin
    143193version = 1
    144194
    145195[domain_root]
    146196pixel_size = 15.0
    147 origin_x = 0.0
    148 origin_y = 0.0
    149 nx = 311
    150 ny = 257
     197origin_x = 0
     198origin_y = 0
     199nx = 3119
     200ny = 2573
    151201buildings_3d = False
    152202dz = 15.0
    153 allow_high_vegetation = True
     203allow_high_vegetation = False
    154204generate_vegetation_patches = True
    155 use_palm_z_axis= True
     205use_palm_z_axis= False
    156206interpolate_terrain = False
    157207domain_parent
     208vegetation_on_roofs = False
     209street_trees = True
    158210
    159211[domain_N02]
     
    170222interpolate_terrain = True
    171223domain_parent = root
     224vegetation_on_roofs = False
     225street_trees = True
  • palm/trunk/SCRIPTS/palm_csd

    r3629 r3668  
    2525# -----------------
    2626# $Id$
     27# Various improvements and bugfixes
     28#
     29# 3629 2018-12-13 12:18:54Z maronga
    2730# Added canopy generator calls. Some improvements
    2831#
     
    7275
    7376#  Definition of settings
    74    global settings_filename_out
     77
    7578   global settings_bridge_width
    7679   global settings_lai_roof_intensive
     
    100103   global global_site
    101104   global global_source
    102    global global_version
     105
     106   global  path_out
     107   global filename_out
     108   global version_out
    103109
    104110
     
    158164   global zt_min
    159165
    160    settings_filename_out = "default"
    161166   settings_bridge_width = 3.0
    162167   settings_lai_roof_intensive = 0.5
     
    185190   global_site = ""
    186191   global_source = ""
    187    global_version = 1
    188 
     192
     193   path_out = ""
     194   version_out = 1
     195   filename_out = "default"
     196   
    189197   domain_names = []
    190198   domain_px = []
     
    260268         global_site = config.get(read_tmp, 'site')
    261269         global_source = config.get(read_tmp, 'source')
    262          global_version = float(config.get(read_tmp, 'version'))
     270
    263271
    264272      if ( read_tmp == 'settings' ):
    265          settings_filename_out = config.get(read_tmp, 'filename_out')
    266273         settings_lai_roof_intensive = config.get(read_tmp, 'lai_roof_intensive')
    267274         settings_lai_roof_extensive = config.get(read_tmp, 'lai_roof_extensive')
     
    273280         settings_lai_alpha = float(config.get(read_tmp, 'lai_alpha'))
    274281         settings_lai_beta = float(config.get(read_tmp, 'lai_beta'))
     282
     283      if ( read_tmp == 'output' ):
     284         path_out = config.get(read_tmp, 'path')         
     285         filename_out = config.get(read_tmp, 'file_out')
     286         version_out = float(config.get(read_tmp, 'version'))
    275287
    276288      if ( read_tmp.split("_")[0] == 'domain' ):
     
    322334         input_file_street_crossings.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_street_crossings'))
    323335         input_file_patch_height.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_patch_height'))
    324          input_file_tree_crown_diameter.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_tree_crown_diameter')) 
     336         
     337         tmp = config.get(read_tmp, 'file_tree_crown_diameter')
     338         if tmp is not None:
     339            input_file_tree_crown_diameter.append(config.get(read_tmp, 'path') + "/" + tmp) 
     340         else:
     341            input_file_tree_crown_diameter.append(None)
    325342         input_file_tree_height.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_tree_height')) 
    326343         input_file_tree_trunk_diameter.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_tree_trunk_diameter')) 
     
    433450#  Calculate indices and input files
    434451   ii.append(input_px.index(domain_px[i]))
    435    filename.append(settings_filename_out + "_" + domain_names[i])
     452   filename.append(path_out + "/" + filename_out + "_" + domain_names[i])
    436453   if domain_parent[i] is not None:
    437454      ii_parent.append(domain_names.index(domain_parent[i]))
     
    447464#  Create NetCDF output file and set global attributes
    448465   nc_create_file(filename[i])
    449    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)
     466   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)
    450467
    451468   
     
    468485del zt_all[:]
    469486
     487print( "Shift terrain height by -" + str(zt_min))
    470488for i in range(0,ndomains):
    471489
     
    475493   y = nc_read_from_file_1d(input_file_y[ii[i]], "y", domain_y0[i], domain_y1[i])
    476494
    477    print( "Shift terrain height by -" + str(zt_min))
    478495   zt = zt - zt_min
     496   
     497   nc_write_global_attribute(filename[i], 'origin_z', float(zt_min))
    479498
    480499#  If necessary, interpolate parent domain terrain height on child domain grid and blend the two
     
    494513
    495514#     Interpolate array and bring to PALM grid of child domain
     515
    496516      zt_ip = interpolate_2d(zt_parent,tmp_x,tmp_y,x,y)
    497517      zt_ip = bring_to_palm_grid(zt_ip,x,y,domain_dz[ii_parent[i]])
     
    580600 
    581601   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]) 
    582    building_type[building_type == 255] = fillvalues["building_type"]
     602   building_type[building_type >= 254] = fillvalues["building_type"]
    583603   building_type = np.where(building_type < 1,defaultvalues["building_type"],building_type)
    584604 
     
    594614      building_type = np.where((building_type == fillvalues["building_type"]) & (buildings_2d != fillvalues["buildings_2d"]),defaultvalues["building_type"],building_type)
    595615      print("Data check #2 " + str(check_arrays_2(buildings_2d,building_type,fillvalues["buildings_2d"],fillvalues["building_type"])))
     616 
    596617 
    597618   nc_write_to_file_2d(filename[i], 'buildings_2d', buildings_2d, datatypes["buildings_2d"],'y','x',fillvalues["buildings_2d"])
     
    606627   nc_write_attribute(filename[i], 'building_id', 'long_name', 'building id')
    607628   nc_write_attribute(filename[i], 'building_id', 'units', '')
    608    nc_write_attribute(filename[i], 'building_id', 'res_orig', domain_px[i])
     629   nc_write_attribute(filename[i], 'building_id', 'res _orig', domain_px[i])
    609630   nc_write_attribute(filename[i], 'building_id', 'coordinates', 'E_UTM N_UTM lon lat')
    610631   nc_write_attribute(filename[i], 'building_id', 'grid_mapping', 'E_UTM N_UTM lon lat')   
     
    616637   nc_write_attribute(filename[i], 'building_type', 'coordinates', 'E_UTM N_UTM lon lat')
    617638   nc_write_attribute(filename[i], 'building_type', 'grid_mapping', 'E_UTM N_UTM lon lat')   
     639 
     640   print("out_1: " + str(np.any(building_type == -2)))
    618641 
    619642del buildings_2d
     
    705728#  #6 Remove water for pixels with buildings
    706729   water_type = np.where((water_type != fillvalues["water_type"]) & (building_type != fillvalues["building_type"]),fillvalues["water_type"],water_type)
    707    
    708 
    709 #  #7 to be removed: set default soil type everywhere
    710    soil_type = np.where((vegetation_type != fillvalues["vegetation_type"]) | (pavement_type != fillvalues["pavement_type"]),defaultvalues["soil_type"],fillvalues["soil_type"])
    711  
     730 
     731   
     732#  Correct vegetation_type when a vegetation height is available and is indicative of low vegeetation
     733   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])   
     734         
     735   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)
     736   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)
    712737 
    713738#  Check for consistency and fill empty fields with default vegetation type
     
    717742      vegetation_type = np.where(consistency_array == 0,defaultvalues["vegetation_type"],vegetation_type)
    718743      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"])
    719    
     744
     745#  #7 to be removed: set default soil type everywhere
     746   soil_type = np.where((vegetation_type != fillvalues["vegetation_type"]) | (pavement_type != fillvalues["pavement_type"]),defaultvalues["soil_type"],fillvalues["soil_type"])
     747       
     748   
     749#  Check for consistency and fill empty fields with default vegetation type
     750   consistency_array, test = check_consistency_3(vegetation_type,pavement_type,soil_type,fillvalues["vegetation_type"],fillvalues["pavement_type"],fillvalues["soil_type"])
     751       
    720752#  Create surface_fraction array
    721753   x = nc_read_from_file_2d_all(filename[i], 'x')
     
    734766   nc_write_attribute(filename[i], 'surface_fraction', 'res_orig', domain_px[i])     
    735767   del surface_fraction   
    736    
    737  
    738 #  Correct vegetation_type when a vegetation height is available and is indicative of low vegeetation
    739    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])   
    740          
    741    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)
    742    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)
    743 
    744768
    745769   nc_write_to_file_2d(filename[i], 'vegetation_type', vegetation_type, datatypes["vegetation_type"],'y','x',fillvalues["vegetation_type"])     
     
    779803
    780804#  pixels with bridges get building_type = 7 = bridge. This does not change the _type setting for the under-bridge
    781 #  area
     805#  area NOTE: when bridges are present the consistency check will fail at the moment
    782806   if domain_3d[i]:
    783807      if np.any(building_type != fillvalues["building_type"]):
     
    853877# Read tree data and create LAD and BAD arrays using the canopy generator
    854878for i in range(0,ndomains): 
     879   lai = nc_read_from_file_2d(input_file_lai[ii[i]], 'Band1', domain_x0[i], domain_x1[i], domain_y0[i], domain_y1[i])   
     880     
     881   vegetation_type = nc_read_from_file_2d_all(filename[i], 'vegetation_type')
     882     
     883   lai = np.where(vegetation_type == fillvalues["vegetation_type"],fillvalues["vegetation_pars"],lai)
     884
     885     
     886   x = nc_read_from_file_2d_all(filename[i], 'x')
     887   y = nc_read_from_file_2d_all(filename[i], 'y') 
     888   nvegetation_pars = np.arange(0,12)
     889   vegetation_pars = np.ones((len(nvegetation_pars),len(y),len(x)))
     890   vegetation_pars[:,:,:] = fillvalues["vegetation_pars"] 
     891     
     892   vegetation_pars[1,:,:] = lai
     893     
     894
     895#  Write out first version of LAI. Will later be overwritten.
     896   nc_write_dimension(filename[i], 'nvegetation_pars', nvegetation_pars, datatypes["nvegetation_pars"])
     897   nc_write_to_file_3d(filename[i], 'vegetation_pars', vegetation_pars, datatypes["vegetation_pars"],'nvegetation_pars','y','x',fillvalues["vegetation_pars"])
     898   nc_write_attribute(filename[i], 'vegetation_pars', 'long_name', 'vegetation_pars')
     899   nc_write_attribute(filename[i], 'vegetation_pars', 'units', '')
     900   nc_write_attribute(filename[i], 'vegetation_pars', 'res_orig', domain_px[i]) 
     901   nc_write_attribute(filename[i], 'vegetation_pars', 'coordinates', 'E_UTM N_UTM lon lat')
     902   nc_write_attribute(filename[i], 'vegetation_pars', 'grid_mapping', 'E_UTM N_UTM lon lat')   
     903
     904   del lai, vegetation_pars, vegetation_type
     905
     906# Read tree data and create LAD and BAD arrays using the canopy generator
     907for i in range(0,ndomains): 
    855908   if domain_street_trees[i]:
    856       lai = nc_read_from_file_2d(input_file_lai[ii[i]], 'Band1', domain_x0[i], domain_x1[i], domain_y0[i], domain_y1[i])   
    857      
    858      
     909
     910      vegetation_pars =  nc_read_from_file_2d_all(filename[i], 'vegetation_pars')
     911
     912      lai = np.copy(vegetation_pars[1,:,:])
     913
     914      x = nc_read_from_file_2d_all(filename[i], 'x')
     915      y = nc_read_from_file_2d_all(filename[i], 'y') 
     916
    859917#     Save lai data as default for low and high vegetation
    860918      lai_low = lai
    861919      lai_high = lai
    862      
    863       x = nc_read_from_file_2d_all(filename[i], 'x')
    864       y = nc_read_from_file_2d_all(filename[i], 'y') 
    865       nvegetation_pars = np.arange(0,11)
    866       vegetation_pars = np.ones((len(nvegetation_pars),len(y),len(x)))
    867       vegetation_pars[:,:,:] = fillvalues["vegetation_pars"] 
    868      
    869       vegetation_pars[1,:,:] = lai
    870      
    871 
    872 #     Write out first version of LAI. Will later be overwritten.
    873       nc_write_dimension(filename[i], 'nvegetation_pars', nvegetation_pars, datatypes["nvegetation_pars"])
    874       nc_write_to_file_3d(filename[i], 'vegetation_pars', vegetation_pars, datatypes["vegetation_pars"],'nvegetation_pars','y','x',fillvalues["vegetation_pars"])
    875       nc_write_attribute(filename[i], 'vegetation_pars', 'long_name', 'vegetation_pars')
    876       nc_write_attribute(filename[i], 'vegetation_pars', 'units', '')
    877       nc_write_attribute(filename[i], 'vegetation_pars', 'res_orig', domain_px[i]) 
    878       nc_write_attribute(filename[i], 'vegetation_pars', 'coordinates', 'E_UTM N_UTM lon lat')
    879       nc_write_attribute(filename[i], 'vegetation_pars', 'grid_mapping', 'E_UTM N_UTM lon lat')   
    880 
    881 
     920   
    882921#     Read all tree parameters from file
    883922      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])   
    884       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])   
     923
     924      if (input_file_tree_crown_diameter[ii[i]] is not None):
     925         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])   
     926         tree_crown_diameter = np.where( (tree_crown_diameter == 0.0) | (tree_crown_diameter == -1.0) ,fillvalues["tree_data"],tree_crown_diameter)
     927      else:
     928         tree_crown_diameter = np.ones((len(y),len(x)))
     929         tree_crown_diameter[:,:] = fillvalues["tree_data"]
     930
     931
    885932      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])         
    886933      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])
     
    889936#     Remove missing values from the data. Reasonable values will be set by the tree generator   
    890937      tree_height = np.where( (tree_height == 0.0) | (tree_height == -1.0) ,fillvalues["tree_data"],tree_height)
    891       tree_crown_diameter = np.where( (tree_crown_diameter == 0.0) | (tree_crown_diameter == -1.0) ,fillvalues["tree_data"],tree_crown_diameter)
    892938      tree_trunk_diameter = np.where( (tree_trunk_diameter == 0.0) | (tree_trunk_diameter == -1.0) ,fillvalues["tree_data"],tree_trunk_diameter)
    893939      tree_type = np.where( (tree_type == 0.0) | (tree_type == -1.0) ,fillvalues["tree_data"],tree_type)
     
    907953#     For zero-height patches, set vegetation_type to short grass and remove these pixels from the patch height array
    908954      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)
    909       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)      
    910       nc_overwrite_to_file_2d(filename[i], 'vegetation_type', vegetation_type) 
     955      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)   
     956     
    911957
    912958      max_tree_height = max(tree_height.flatten())
     
    914960     
    915961      if ( (max_tree_height != fillvalues["tree_data"]) | (max_patch_height == fillvalues["tree_data"]) ):
    916            
     962
    917963         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"])
    918964 
     
    925971            bad[k,:,:] = np.where(buildings_2d == fillvalues["buildings_2d"],bad[k,:,:],fillvalues["tree_data"])
    926972            tree_ids[k,:,:] = np.where(buildings_2d == fillvalues["buildings_2d"],tree_ids[k,:,:],fillvalues["tree_data"])     
    927  
     973
     974         del buildings_2d
    928975 
    929976         nc_write_dimension(filename[i], 'zlad', zlad, datatypes["tree_data"])
     
    949996         nc_write_attribute(filename[i], 'tree_id', 'grid_mapping', 'E_UTM N_UTM lon lat')           
    950997           
    951          del buildings_2d, lai, lad, bad, tree_ids, zlad
    952 
    953       del vegetation_pars, tree_height, tree_crown_diameter, tree_trunk_diameter, tree_type, patch_height, vegetation_type, x, y
     998         del lai, lad, bad, tree_ids, zlad
     999
     1000      del vegetation_pars, tree_height, tree_crown_diameter, tree_trunk_diameter, tree_type, patch_height, x, y
    9541001
    9551002
     
    9611008      vegetation_type = nc_read_from_file_2d_all(filename[i], 'vegetation_type')
    9621009      lad = nc_read_from_file_3d_all(filename[i], 'lad')
     1010      zlad = nc_read_from_file_1d_all(filename[i], 'zlad')
    9631011      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])
    9641012      vegetation_pars = nc_read_from_file_3d_all(filename[i], 'vegetation_pars')
    9651013      lai = vegetation_pars[1,:,:]
    9661014
    967 #     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
    968 #     the lai_high is reset to a low vegetation value
    969       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"])
    970 
    971 #     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?!
    972       lai_high = np.where( (lai_high != fillvalues["tree_data"]) & (lai == fillvalues["tree_data"]), settings_lai_high_default, lai_high)
    973 
    974 #     Define a patch height whereever it is missing, but where a high vegetation LAI was set
     1015
     1016#     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
     1017      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"])
     1018
     1019#     Now, assign either the default LAI for high vegetation or the value stored in the LAI array.
     1020      lai_high = np.where( (lai_high != fillvalues["tree_data"]) & (lai == fillvalues["tree_data"]), settings_lai_high_default, lai)
     1021
     1022#     Define a patch height wherever it is missing, but where a high vegetation LAI was set
    9751023      patch_height = np.where ( (lai_high != fillvalues["tree_data"]) & (patch_height == fillvalues["tree_data"]), settings_patch_height_default, patch_height)
    9761024
     
    9781026      patch_height = np.where ( (lad[0,:,:] != fillvalues["tree_data"]), fillvalues["tree_data"], patch_height)
    9791027
     1028#     Remove patch heights that have no lai_high value
     1029      patch_height = np.where ( (lai_high == fillvalues["tree_data"]), fillvalues["tree_data"], patch_height)
     1030
    9801031#     For missing LAI values, set either the high vegetation default or the low vegetation default   
    981       lai_high = np.where( (patch_height > 2.0) & (lai_high == fillvalues["tree_data"]),settings_lai_high_default,lai_high)
    982       lai_high = np.where( (patch_height <= 2.0) & (lai_high == fillvalues["tree_data"]),settings_lai_low_default,lai_high)
    983    
    984    
     1032      lai_high = np.where( (patch_height > 2.0) & (patch_height != fillvalues["tree_data"]) & (lai_high == fillvalues["tree_data"]),settings_lai_high_default,lai_high)
     1033      lai_high = np.where( (patch_height <= 2.0) & (patch_height != fillvalues["tree_data"]) & (lai_high == fillvalues["tree_data"]),settings_lai_low_default,lai_high)
     1034         
    9851035      if ( max(patch_height.flatten()) >= (2.0 * domain_dz[i]) ):
    9861036         print("    start calculating LAD (this might take some time)")
    987          lad_patch, patch_nz, status = process_patch(domain_dz[i],patch_height,lai_high,settings_lai_alpha,settings_lai_beta)
     1037         lad_patch, patch_nz, status = process_patch(domain_dz[i],patch_height,max(zlad),lai_high,settings_lai_alpha,settings_lai_beta)
    9881038
    9891039         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,:,:] )
     
    9911041#     Remove high vegetation wherever it is replaced by a leaf area density. This should effectively remove all high vegetation pixels
    9921042      vegetation_type = np.where((lad[0,:,:] != fillvalues["tree_data"]) & (vegetation_type != fillvalues["vegetation_type"]),3,vegetation_type)
    993  
    994  
     1043             
    9951044#     If desired, remove all high vegetation. TODO: check if this is still necessary
    9961045      if domain_high_vegetation[i]:
     
    9981047
    9991048 
    1000 #     Remove LAI from pixels where a LAD was set
    1001       lai_low = np.where( (lad[0,:,:] == fillvalues["tree_data"]), lai,fillvalues["tree_data"])
     1049#     Set default low LAI for pixels with an LAD (short grass below trees)
     1050      lai_low = np.where( (lad[0,:,:] == fillvalues["tree_data"]), lai, settings_lai_low_default)
    10021051   
    1003 #     Fill low vegetation pixels without LAI set with default value
    1004       lai_low = np.where( (lai == fillvalues["tree_data"]) & (vegetation_type != fillvalues["vegetation_type"] ), settings_lai_low_default, lai)
    1005 
    1006    
     1052#     Fill low vegetation pixels without LAI set or with LAI = 0 with default value
     1053      lai_low = np.where( ( (lai_low == fillvalues["tree_data"]) | (lai_low == 0.0) ) & (vegetation_type != fillvalues["vegetation_type"] ), settings_lai_low_default, lai_low)
     1054
     1055#     Remove lai for pixels that have no vegetation_type
     1056      lai_low = np.where( vegetation_type != fillvalues["vegetation_type"], lai_low, fillvalues["tree_data"])
     1057       
    10071058#     Overwrite lai in vegetation_parameters 
    1008 
    1009       vegetation_pars[1,:,:] = lai_low
     1059      vegetation_pars[1,:,:] = np.copy(lai_low)
    10101060      nc_overwrite_to_file_3d(filename[i], 'vegetation_pars', vegetation_pars)
    10111061
    10121062#     Overwrite lad array
    10131063      nc_overwrite_to_file_3d(filename[i], 'lad', lad)
    1014      
    1015 #     Overwrite vegetation_type
    1016       nc_overwrite_to_file_2d(filename[i], 'vegetation_type', vegetation_type)   
    1017      
    1018      
    1019       del vegetation_type, lad, lai, patch_height, vegetation_pars
     1064           
     1065      nc_overwrite_to_file_2d(filename[i], 'vegetation_type', vegetation_type) 
     1066
     1067     
     1068      del vegetation_type, lad, lai, patch_height, vegetation_pars, zlad
     1069
     1070# Final consistency check
     1071for i in range(0,ndomains):
     1072   vegetation_type = nc_read_from_file_2d_all(filename[i], 'vegetation_type') 
     1073   pavement_type = nc_read_from_file_2d_all(filename[i], 'pavement_type') 
     1074   building_type = nc_read_from_file_2d_all(filename[i], 'building_type') 
     1075   water_type = nc_read_from_file_2d_all(filename[i], 'water_type')     
     1076   soil_type = nc_read_from_file_2d_all(filename[i], 'soil_type') 
     1077
     1078#  Check for consistency and fill empty fields with default vegetation type
     1079   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"])
     1080
     1081#  Check for consistency and fill empty fields with default vegetation type
     1082   consistency_array, test = check_consistency_3(vegetation_type,pavement_type,soil_type,fillvalues["vegetation_type"],fillvalues["pavement_type"],fillvalues["soil_type"])
     1083
     1084   surface_fraction = nc_read_from_file_3d_all(filename[i], 'surface_fraction')       
     1085   surface_fraction[0,:,:] = np.where(vegetation_type != fillvalues["vegetation_type"], 1.0, 0.0)
     1086   surface_fraction[1,:,:] = np.where(pavement_type != fillvalues["pavement_type"], 1.0, 0.0)
     1087   surface_fraction[2,:,:] = np.where(water_type != fillvalues["water_type"], 1.0, 0.0)
     1088   nc_overwrite_to_file_3d(filename[i], 'surface_fraction', surface_fraction)
     1089       
     1090   del vegetation_type, pavement_type, building_type, water_type, soil_type
  • palm/trunk/SCRIPTS/palm_csd_files/palm_csd_canopy_generator.py

    r3629 r3668  
    2525# -----------------
    2626# $Id$
     27# Various improvements and bugfixes
     28#
     29# 3629 2018-12-13 12:18:54Z maronga
    2730# Initial revision
    2831#
     
    7275            if valid_pixels[j,i]:
    7376               tree_id_counter = tree_id_counter + 1
    74                print("   Processing tree No " +  str(tree_id_counter) + " ...", end="")
    75 
     77#               print("   Processing tree No " +  str(tree_id_counter) + " ...", end="")
    7678               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)
    7779               if ( np.any(lad_loc) != fill ):
     
    101103
    102104
    103                   if ( status == 0 ):
    104                      status_char = " ok."
    105                   else:
    106                      status_char = " skipped."
    107                   print(status_char)
     105#                  if ( status == 0 ):
     106#                     status_char = " ok."
     107#                  else:
     108#                     status_char = " skipped."
     109#                  print(status_char)
    108110
    109111               del lad_loc, x_loc, y_loc, z_loc, status
     
    131133   default_trees = []
    132134                                                #1   #2   #3   #4    #5   #6   #7   #8     #9
    133    default_trees.append(tree("Default",         1.0, 1.0, 4.0, 12.0, 2.0, 0.1, 0.6, 0.025, 0.35))
    134    default_trees.append(tree("Abies",           3.0, 1.0, 4.0, 12.0, 2.0, 0.1, 0.6, 0.025, 0.80))
    135    default_trees.append(tree("Acer",            1.0, 1.0, 7.0, 12.0, 2.0, 0.1, 0.6, 0.025, 0.80))
    136    default_trees.append(tree("Aesculus",        1.0, 1.0, 7.0, 12.0, 2.0, 0.1, 0.6, 0.025, 1.00))
    137    default_trees.append(tree("Ailanthus",       1.0, 1.0, 8.5, 13.5, 2.0, 0.1, 0.6, 0.025, 1.30))
    138    default_trees.append(tree("Alnus",           3.0, 1.0, 6.0, 16.0, 2.0, 0.1, 0.6, 0.025, 1.20))
    139    default_trees.append(tree("Amelanchier",     1.0, 1.0, 3.0,  4.0, 2.0, 0.1, 0.6, 0.025, 1.20))
    140    default_trees.append(tree("Betula",          1.0, 1.0, 6.0, 14.0, 2.0, 0.1, 0.6, 0.025, 0.30))
    141    default_trees.append(tree("Buxus",           1.0, 1.0, 4.0,  4.0, 2.0, 0.1, 0.6, 0.025, 0.90))
    142    default_trees.append(tree("Calocedrus",      3.0, 1.0, 5.0, 10.0, 2.0, 0.1, 0.6, 0.025, 0.50))
    143    default_trees.append(tree("Caragana",        1.0, 1.0, 3.5,  6.0, 2.0, 0.1, 0.6, 0.025, 0.90))
    144    default_trees.append(tree("Carpinus",        1.0, 1.0, 6.0, 10.0, 2.0, 0.1, 0.6, 0.025, 0.70))
    145    default_trees.append(tree("Carya",           1.0, 1.0, 5.0, 17.0, 2.0, 0.1, 0.6, 0.025, 0.80))
    146    default_trees.append(tree("Castanea",        1.0, 1.0, 4.5,  7.0, 2.0, 0.1, 0.6, 0.025, 0.80))
    147    default_trees.append(tree("Catalpa",         1.0, 1.0, 5.5,  6.5, 2.0, 0.1, 0.6, 0.025, 0.70))
    148    default_trees.append(tree("Cedrus",          1.0, 1.0, 8.0, 13.0, 2.0, 0.1, 0.6, 0.025, 0.80))
    149    default_trees.append(tree("Celtis",          1.0, 1.0, 6.0,  9.0, 2.0, 0.1, 0.6, 0.025, 0.80))
    150    default_trees.append(tree("Cercidiphyllum",  1.0, 1.0, 3.0,  6.5, 2.0, 0.1, 0.6, 0.025, 0.80))
    151    default_trees.append(tree("Cercis",          1.0, 1.0, 2.5,  7.5, 2.0, 0.1, 0.6, 0.025, 0.90))
    152    default_trees.append(tree("Chamaecyparis",   5.0, 1.0, 3.5,  9.0, 2.0, 0.1, 0.6, 0.025, 0.70))
    153    default_trees.append(tree("Cladrastis",      1.0, 1.0, 5.0, 10.0, 2.0, 0.1, 0.6, 0.025, 0.80))
    154    default_trees.append(tree("Cornus",          1.0, 1.0, 4.5,  6.5, 2.0, 0.1, 0.6, 0.025, 1.20))
    155    default_trees.append(tree("Corylus",         1.0, 1.0, 5.0,  9.0, 2.0, 0.1, 0.6, 0.025, 0.40))
    156    default_trees.append(tree("Cotinus",         1.0, 1.0, 4.0,  4.0, 2.0, 0.1, 0.6, 0.025, 0.70))
    157    default_trees.append(tree("Crataegus",       3.0, 1.0, 3.5,  6.0, 2.0, 0.1, 0.6, 0.025, 1.40))
    158    default_trees.append(tree("Cryptomeria",     3.0, 1.0, 5.0, 10.0, 2.0, 0.1, 0.6, 0.025, 0.50))
    159    default_trees.append(tree("Cupressocyparis", 3.0, 1.0, 3.0,  8.0, 2.0, 0.1, 0.6, 0.025, 0.40))
    160    default_trees.append(tree("Cupressus",       3.0, 1.0, 5.0,  7.0, 2.0, 0.1, 0.6, 0.025, 0.40))
    161    default_trees.append(tree("Cydonia",         1.0, 1.0, 2.0,  3.0, 2.0, 0.1, 0.6, 0.025, 0.90))
    162    default_trees.append(tree("Davidia",         1.0, 1.0,10.0, 14.0, 2.0, 0.1, 0.6, 0.025, 0.40))
    163    default_trees.append(tree("Elaeagnus",       1.0, 1.0, 6.5,  6.0, 2.0, 0.1, 0.6, 0.025, 1.20))
    164    default_trees.append(tree("Euodia",          1.0, 1.0, 4.5,  6.0, 2.0, 0.1, 0.6, 0.025, 0.90))
    165    default_trees.append(tree("Euonymus",        1.0, 1.0, 4.5,  6.0, 2.0, 0.1, 0.6, 0.025, 0.60))
    166    default_trees.append(tree("Fagus",           1.0, 1.0,10.0, 12.5, 2.0, 0.1, 0.6, 0.025, 0.50))
    167    default_trees.append(tree("Fraxinus",        1.0, 1.0, 5.5, 10.5, 2.0, 0.1, 0.6, 0.025, 1.60))
    168    default_trees.append(tree("Ginkgo",          3.0, 1.0, 4.0,  8.5, 2.0, 0.1, 0.6, 0.025, 0.80))
    169    default_trees.append(tree("Gleditsia",       1.0, 1.0, 6.5, 10.5, 2.0, 0.1, 0.6, 0.025, 0.60))
    170    default_trees.append(tree("Gymnocladus",     1.0, 1.0, 5.5, 10.0, 2.0, 0.1, 0.6, 0.025, 0.80))
    171    default_trees.append(tree("Hippophae",       1.0, 1.0, 9.5,  8.5, 2.0, 0.1, 0.6, 0.025, 0.80))
    172    default_trees.append(tree("Ilex",            1.0, 1.0, 4.0,  7.5, 2.0, 0.1, 0.6, 0.025, 0.80))
    173    default_trees.append(tree("Juglans",         1.0, 1.0, 7.0,  9.0, 2.0, 0.1, 0.6, 0.025, 0.50))
    174    default_trees.append(tree("Juniperus",       5.0, 1.0, 3.0,  7.0, 2.0, 0.1, 0.6, 0.025, 0.90))
    175    default_trees.append(tree("Koelreuteria",    1.0, 1.0, 3.5,  5.5, 2.0, 0.1, 0.6, 0.025, 0.50))
    176    default_trees.append(tree("Laburnum",        1.0, 1.0, 3.0,  6.0, 2.0, 0.1, 0.6, 0.025, 0.60))
    177    default_trees.append(tree("Larix",           3.0, 1.0, 7.0, 16.5, 2.0, 0.1, 0.6, 0.025, 0.60))
    178    default_trees.append(tree("Ligustrum",       1.0, 1.0, 3.0,  6.0, 2.0, 0.1, 0.6, 0.025, 1.10))
    179    default_trees.append(tree("Liquidambar",     3.0, 1.0, 3.0,  7.0, 2.0, 0.1, 0.6, 0.025, 0.30))
    180    default_trees.append(tree("Liriodendron",    3.0, 1.0, 4.5,  9.5, 2.0, 0.1, 0.6, 0.025, 0.50))
    181    default_trees.append(tree("Lonicera",        1.0, 1.0, 7.0,  9.0, 2.0, 0.1, 0.6, 0.025, 0.70))
    182    default_trees.append(tree("Magnolia",        1.0, 1.0, 3.0,  5.0, 2.0, 0.1, 0.6, 0.025, 0.60))
    183    default_trees.append(tree("Malus",           1.0, 1.0, 4.5,  5.0, 2.0, 0.1, 0.6, 0.025, 0.30))
    184    default_trees.append(tree("Metasequoia",     5.0, 1.0, 4.5, 12.0, 2.0, 0.1, 0.6, 0.025, 0.50))
    185    default_trees.append(tree("Morus",           1.0, 1.0, 7.5, 11.5, 2.0, 0.1, 0.6, 0.025, 1.00))
    186    default_trees.append(tree("Ostrya",          1.0, 1.0, 2.0,  6.0, 2.0, 0.1, 0.6, 0.025, 1.00))
    187    default_trees.append(tree("Parrotia",        1.0, 1.0, 7.0,  7.0, 2.0, 0.1, 0.6, 0.025, 0.30))
    188    default_trees.append(tree("Paulownia",       1.0, 1.0, 4.0,  8.0, 2.0, 0.1, 0.6, 0.025, 0.40))
    189    default_trees.append(tree("Phellodendron",   1.0, 1.0,13.5, 13.5, 2.0, 0.1, 0.6, 0.025, 0.50))
    190    default_trees.append(tree("Picea",           3.0, 1.0, 3.0, 13.0, 2.0, 0.1, 0.6, 0.025, 0.90))
    191    default_trees.append(tree("Pinus",           3.0, 1.0, 6.0, 16.0, 2.0, 0.1, 0.6, 0.025, 0.80))
    192    default_trees.append(tree("Platanus",        1.0, 1.0,10.0, 14.5, 2.0, 0.1, 0.6, 0.025, 1.10))
    193    default_trees.append(tree("Populus",         1.0, 1.0, 9.0, 20.0, 2.0, 0.1, 0.6, 0.025, 1.40))
    194    default_trees.append(tree("Prunus",          1.0, 1.0, 5.0,  7.0, 2.0, 0.1, 0.6, 0.025, 1.60))
    195    default_trees.append(tree("Pseudotsuga",     3.0, 1.0, 6.0, 17.5, 2.0, 0.1, 0.6, 0.025, 0.70))
    196    default_trees.append(tree("Ptelea",          1.0, 1.0, 5.0,  4.0, 2.0, 0.1, 0.6, 0.025, 1.10))
    197    default_trees.append(tree("Pterocaria",      1.0, 1.0,10.0, 12.0, 2.0, 0.1, 0.6, 0.025, 0.50))
    198    default_trees.append(tree("Pterocarya",      1.0, 1.0,11.5, 14.5, 2.0, 0.1, 0.6, 0.025, 1.60))
    199    default_trees.append(tree("Pyrus",           3.0, 1.0, 3.0,  6.0, 2.0, 0.1, 0.6, 0.025, 1.80))
     135   default_trees.append(tree("Default",         1.0, 1.0, 4.0, 12.0, 3.0, 0.8, 0.6, 0.025, 0.35))
     136   default_trees.append(tree("Abies",           3.0, 1.0, 4.0, 12.0, 3.0, 0.8, 0.6, 0.025, 0.80))
     137   default_trees.append(tree("Acer",            1.0, 1.0, 7.0, 12.0, 3.0, 0.8, 0.6, 0.025, 0.80))
     138   default_trees.append(tree("Aesculus",        1.0, 1.0, 7.0, 12.0, 3.0, 0.8, 0.6, 0.025, 1.00))
     139   default_trees.append(tree("Ailanthus",       1.0, 1.0, 8.5, 13.5, 3.0, 0.8, 0.6, 0.025, 1.30))
     140   default_trees.append(tree("Alnus",           3.0, 1.0, 6.0, 16.0, 3.0, 0.8, 0.6, 0.025, 1.20))
     141   default_trees.append(tree("Amelanchier",     1.0, 1.0, 3.0,  4.0, 3.0, 0.8, 0.6, 0.025, 1.20))
     142   default_trees.append(tree("Betula",          1.0, 1.0, 6.0, 14.0, 3.0, 0.8, 0.6, 0.025, 0.30))
     143   default_trees.append(tree("Buxus",           1.0, 1.0, 4.0,  4.0, 3.0, 0.8, 0.6, 0.025, 0.90))
     144   default_trees.append(tree("Calocedrus",      3.0, 1.0, 5.0, 10.0, 3.0, 0.8, 0.6, 0.025, 0.50))
     145   default_trees.append(tree("Caragana",        1.0, 1.0, 3.5,  6.0, 3.0, 0.8, 0.6, 0.025, 0.90))
     146   default_trees.append(tree("Carpinus",        1.0, 1.0, 6.0, 10.0, 3.0, 0.8, 0.6, 0.025, 0.70))
     147   default_trees.append(tree("Carya",           1.0, 1.0, 5.0, 17.0, 3.0, 0.8, 0.6, 0.025, 0.80))
     148   default_trees.append(tree("Castanea",        1.0, 1.0, 4.5,  7.0, 3.0, 0.8, 0.6, 0.025, 0.80))
     149   default_trees.append(tree("Catalpa",         1.0, 1.0, 5.5,  6.5, 3.0, 0.8, 0.6, 0.025, 0.70))
     150   default_trees.append(tree("Cedrus",          1.0, 1.0, 8.0, 13.0, 3.0, 0.8, 0.6, 0.025, 0.80))
     151   default_trees.append(tree("Celtis",          1.0, 1.0, 6.0,  9.0, 3.0, 0.8, 0.6, 0.025, 0.80))
     152   default_trees.append(tree("Cercidiphyllum",  1.0, 1.0, 3.0,  6.5, 3.0, 0.8, 0.6, 0.025, 0.80))
     153   default_trees.append(tree("Cercis",          1.0, 1.0, 2.5,  7.5, 3.0, 0.8, 0.6, 0.025, 0.90))
     154   default_trees.append(tree("Chamaecyparis",   5.0, 1.0, 3.5,  9.0, 3.0, 0.8, 0.6, 0.025, 0.70))
     155   default_trees.append(tree("Cladrastis",      1.0, 1.0, 5.0, 10.0, 3.0, 0.8, 0.6, 0.025, 0.80))
     156   default_trees.append(tree("Cornus",          1.0, 1.0, 4.5,  6.5, 3.0, 0.8, 0.6, 0.025, 1.20))
     157   default_trees.append(tree("Corylus",         1.0, 1.0, 5.0,  9.0, 3.0, 0.8, 0.6, 0.025, 0.40))
     158   default_trees.append(tree("Cotinus",         1.0, 1.0, 4.0,  4.0, 3.0, 0.8, 0.6, 0.025, 0.70))
     159   default_trees.append(tree("Crataegus",       3.0, 1.0, 3.5,  6.0, 3.0, 0.8, 0.6, 0.025, 1.40))
     160   default_trees.append(tree("Cryptomeria",     3.0, 1.0, 5.0, 10.0, 3.0, 0.8, 0.6, 0.025, 0.50))
     161   default_trees.append(tree("Cupressocyparis", 3.0, 1.0, 3.0,  8.0, 3.0, 0.8, 0.6, 0.025, 0.40))
     162   default_trees.append(tree("Cupressus",       3.0, 1.0, 5.0,  7.0, 3.0, 0.8, 0.6, 0.025, 0.40))
     163   default_trees.append(tree("Cydonia",         1.0, 1.0, 2.0,  3.0, 3.0, 0.8, 0.6, 0.025, 0.90))
     164   default_trees.append(tree("Davidia",         1.0, 1.0,10.0, 14.0, 3.0, 0.8, 0.6, 0.025, 0.40))
     165   default_trees.append(tree("Elaeagnus",       1.0, 1.0, 6.5,  6.0, 3.0, 0.8, 0.6, 0.025, 1.20))
     166   default_trees.append(tree("Euodia",          1.0, 1.0, 4.5,  6.0, 3.0, 0.8, 0.6, 0.025, 0.90))
     167   default_trees.append(tree("Euonymus",        1.0, 1.0, 4.5,  6.0, 3.0, 0.8, 0.6, 0.025, 0.60))
     168   default_trees.append(tree("Fagus",           1.0, 1.0,10.0, 12.5, 3.0, 0.8, 0.6, 0.025, 0.50))
     169   default_trees.append(tree("Fraxinus",        1.0, 1.0, 5.5, 10.5, 3.0, 0.8, 0.6, 0.025, 1.60))
     170   default_trees.append(tree("Ginkgo",          3.0, 1.0, 4.0,  8.5, 3.0, 0.8, 0.6, 0.025, 0.80))
     171   default_trees.append(tree("Gleditsia",       1.0, 1.0, 6.5, 10.5, 3.0, 0.8, 0.6, 0.025, 0.60))
     172   default_trees.append(tree("Gymnocladus",     1.0, 1.0, 5.5, 10.0, 3.0, 0.8, 0.6, 0.025, 0.80))
     173   default_trees.append(tree("Hippophae",       1.0, 1.0, 9.5,  8.5, 3.0, 0.8, 0.6, 0.025, 0.80))
     174   default_trees.append(tree("Ilex",            1.0, 1.0, 4.0,  7.5, 3.0, 0.8, 0.6, 0.025, 0.80))
     175   default_trees.append(tree("Juglans",         1.0, 1.0, 7.0,  9.0, 3.0, 0.8, 0.6, 0.025, 0.50))
     176   default_trees.append(tree("Juniperus",       5.0, 1.0, 3.0,  7.0, 3.0, 0.8, 0.6, 0.025, 0.90))
     177   default_trees.append(tree("Koelreuteria",    1.0, 1.0, 3.5,  5.5, 3.0, 0.8, 0.6, 0.025, 0.50))
     178   default_trees.append(tree("Laburnum",        1.0, 1.0, 3.0,  6.0, 3.0, 0.8, 0.6, 0.025, 0.60))
     179   default_trees.append(tree("Larix",           3.0, 1.0, 7.0, 16.5, 3.0, 0.8, 0.6, 0.025, 0.60))
     180   default_trees.append(tree("Ligustrum",       1.0, 1.0, 3.0,  6.0, 3.0, 0.8, 0.6, 0.025, 1.10))
     181   default_trees.append(tree("Liquidambar",     3.0, 1.0, 3.0,  7.0, 3.0, 0.8, 0.6, 0.025, 0.30))
     182   default_trees.append(tree("Liriodendron",    3.0, 1.0, 4.5,  9.5, 3.0, 0.8, 0.6, 0.025, 0.50))
     183   default_trees.append(tree("Lonicera",        1.0, 1.0, 7.0,  9.0, 3.0, 0.8, 0.6, 0.025, 0.70))
     184   default_trees.append(tree("Magnolia",        1.0, 1.0, 3.0,  5.0, 3.0, 0.8, 0.6, 0.025, 0.60))
     185   default_trees.append(tree("Malus",           1.0, 1.0, 4.5,  5.0, 3.0, 0.8, 0.6, 0.025, 0.30))
     186   default_trees.append(tree("Metasequoia",     5.0, 1.0, 4.5, 12.0, 3.0, 0.8, 0.6, 0.025, 0.50))
     187   default_trees.append(tree("Morus",           1.0, 1.0, 7.5, 11.5, 3.0, 0.8, 0.6, 0.025, 1.00))
     188   default_trees.append(tree("Ostrya",          1.0, 1.0, 2.0,  6.0, 3.0, 0.8, 0.6, 0.025, 1.00))
     189   default_trees.append(tree("Parrotia",        1.0, 1.0, 7.0,  7.0, 3.0, 0.8, 0.6, 0.025, 0.30))
     190   default_trees.append(tree("Paulownia",       1.0, 1.0, 4.0,  8.0, 3.0, 0.8, 0.6, 0.025, 0.40))
     191   default_trees.append(tree("Phellodendron",   1.0, 1.0,13.5, 13.5, 3.0, 0.8, 0.6, 0.025, 0.50))
     192   default_trees.append(tree("Picea",           3.0, 1.0, 3.0, 13.0, 3.0, 0.8, 0.6, 0.025, 0.90))
     193   default_trees.append(tree("Pinus",           3.0, 1.0, 6.0, 16.0, 3.0, 0.8, 0.6, 0.025, 0.80))
     194   default_trees.append(tree("Platanus",        1.0, 1.0,10.0, 14.5, 3.0, 0.8, 0.6, 0.025, 1.10))
     195   default_trees.append(tree("Populus",         1.0, 1.0, 9.0, 20.0, 3.0, 0.8, 0.6, 0.025, 1.40))
     196   default_trees.append(tree("Prunus",          1.0, 1.0, 5.0,  7.0, 3.0, 0.8, 0.6, 0.025, 1.60))
     197   default_trees.append(tree("Pseudotsuga",     3.0, 1.0, 6.0, 17.5, 3.0, 0.8, 0.6, 0.025, 0.70))
     198   default_trees.append(tree("Ptelea",          1.0, 1.0, 5.0,  4.0, 3.0, 0.8, 0.6, 0.025, 1.10))
     199   default_trees.append(tree("Pterocaria",      1.0, 1.0,10.0, 12.0, 3.0, 0.8, 0.6, 0.025, 0.50))
     200   default_trees.append(tree("Pterocarya",      1.0, 1.0,11.5, 14.5, 3.0, 0.8, 0.6, 0.025, 1.60))
     201   default_trees.append(tree("Pyrus",           3.0, 1.0, 3.0,  6.0, 3.0, 0.8, 0.6, 0.025, 1.80))
    200202   default_trees.append(tree("Quercus",         1.0, 1.0, 8.0, 14.0, 3.1, 0.1, 0.6, 0.025, 0.40)) #
    201    default_trees.append(tree("Rhamnus",         1.0, 1.0, 4.5,  4.5, 2.0, 0.1, 0.6, 0.025, 1.30))
    202    default_trees.append(tree("Rhus",            1.0, 1.0, 7.0,  5.5, 2.0, 0.1, 0.6, 0.025, 0.50))
    203    default_trees.append(tree("Robinia",         1.0, 1.0, 4.5, 13.5, 2.0, 0.1, 0.6, 0.025, 0.50))
    204    default_trees.append(tree("Salix",           1.0, 1.0, 7.0, 14.0, 2.0, 0.1, 0.6, 0.025, 1.10))
    205    default_trees.append(tree("Sambucus",        1.0, 1.0, 8.0,  6.0, 2.0, 0.1, 0.6, 0.025, 1.40))
    206    default_trees.append(tree("Sasa",            1.0, 1.0,10.0, 25.0, 2.0, 0.1, 0.6, 0.025, 0.60))
    207    default_trees.append(tree("Sequoiadendron",  5.0, 1.0, 5.5, 10.5, 2.0, 0.1, 0.6, 0.025, 1.60))
    208    default_trees.append(tree("Sophora",         1.0, 1.0, 7.5, 10.0, 2.0, 0.1, 0.6, 0.025, 1.40))
    209    default_trees.append(tree("Sorbus",          1.0, 1.0, 4.0,  7.0, 2.0, 0.1, 0.6, 0.025, 1.10))
    210    default_trees.append(tree("Syringa",         1.0, 1.0, 4.5,  5.0, 2.0, 0.1, 0.6, 0.025, 0.60))
    211    default_trees.append(tree("Tamarix",         1.0, 1.0, 6.0,  7.0, 2.0, 0.1, 0.6, 0.025, 0.50))
    212    default_trees.append(tree("Taxodium",        5.0, 1.0, 6.0, 16.5, 2.0, 0.1, 0.6, 0.025, 0.60))
    213    default_trees.append(tree("Taxus",           2.0, 1.0, 5.0,  7.5, 2.0, 0.1, 0.6, 0.025, 1.50))
    214    default_trees.append(tree("Thuja",           3.0, 1.0, 3.5,  9.0, 2.0, 0.1, 0.6, 0.025, 0.70))
    215    default_trees.append(tree("Tilia",           3.0, 1.0, 7.0, 12.5, 2.0, 0.1, 0.6, 0.025, 0.70))
    216    default_trees.append(tree("Tsuga",           3.0, 1.0, 6.0, 10.5, 2.0, 0.1, 0.6, 0.025, 1.10))
    217    default_trees.append(tree("Ulmus",           1.0, 1.0, 7.5, 14.0, 2.0, 0.1, 0.6, 0.025, 0.80))
    218    default_trees.append(tree("Zelkova",         1.0, 1.0, 4.0,  5.5, 2.0, 0.1, 0.6, 0.025, 1.20))
    219    default_trees.append(tree("Zenobia",         1.0, 1.0, 5.0,  5.0, 2.0, 0.1, 0.6, 0.025, 0.40))
     203   default_trees.append(tree("Rhamnus",         1.0, 1.0, 4.5,  4.5, 3.0, 0.8, 0.6, 0.025, 1.30))
     204   default_trees.append(tree("Rhus",            1.0, 1.0, 7.0,  5.5, 3.0, 0.8, 0.6, 0.025, 0.50))
     205   default_trees.append(tree("Robinia",         1.0, 1.0, 4.5, 13.5, 3.0, 0.8, 0.6, 0.025, 0.50))
     206   default_trees.append(tree("Salix",           1.0, 1.0, 7.0, 14.0, 3.0, 0.8, 0.6, 0.025, 1.10))
     207   default_trees.append(tree("Sambucus",        1.0, 1.0, 8.0,  6.0, 3.0, 0.8, 0.6, 0.025, 1.40))
     208   default_trees.append(tree("Sasa",            1.0, 1.0,10.0, 25.0, 3.0, 0.8, 0.6, 0.025, 0.60))
     209   default_trees.append(tree("Sequoiadendron",  5.0, 1.0, 5.5, 10.5, 3.0, 0.8, 0.6, 0.025, 1.60))
     210   default_trees.append(tree("Sophora",         1.0, 1.0, 7.5, 10.0, 3.0, 0.8, 0.6, 0.025, 1.40))
     211   default_trees.append(tree("Sorbus",          1.0, 1.0, 4.0,  7.0, 3.0, 0.8, 0.6, 0.025, 1.10))
     212   default_trees.append(tree("Syringa",         1.0, 1.0, 4.5,  5.0, 3.0, 0.8, 0.6, 0.025, 0.60))
     213   default_trees.append(tree("Tamarix",         1.0, 1.0, 6.0,  7.0, 3.0, 0.8, 0.6, 0.025, 0.50))
     214   default_trees.append(tree("Taxodium",        5.0, 1.0, 6.0, 16.5, 3.0, 0.8, 0.6, 0.025, 0.60))
     215   default_trees.append(tree("Taxus",           2.0, 1.0, 5.0,  7.5, 3.0, 0.8, 0.6, 0.025, 1.50))
     216   default_trees.append(tree("Thuja",           3.0, 1.0, 3.5,  9.0, 3.0, 0.8, 0.6, 0.025, 0.70))
     217   default_trees.append(tree("Tilia",           3.0, 1.0, 7.0, 12.5, 3.0, 0.8, 0.6, 0.025, 0.70))
     218   default_trees.append(tree("Tsuga",           3.0, 1.0, 6.0, 10.5, 3.0, 0.8, 0.6, 0.025, 1.10))
     219   default_trees.append(tree("Ulmus",           1.0, 1.0, 7.5, 14.0, 3.0, 0.8, 0.6, 0.025, 0.80))
     220   default_trees.append(tree("Zelkova",         1.0, 1.0, 4.0,  5.5, 3.0, 0.8, 0.6, 0.025, 1.20))
     221   default_trees.append(tree("Zenobia",         1.0, 1.0, 5.0,  5.0, 3.0, 0.8, 0.6, 0.025, 0.40))
    220222   
    221223#  Define fill values
     
    237239      if (season == "summer"):
    238240         tree_lai = default_trees[tree_type].lai_summer
     241
    239242      else:
    240243         tree_lai = default_trees[tree_type].lai_winter     
     
    279282   
    280283   lad_max = tree_lai / (lad_max_part_1[0] + lad_max_part_2[0])
    281 
    282284   
    283285   
     
    298300   z = np.arange(0,nz*dz,dz)
    299301   z[1:] = z[1:] - 0.5 * dz
    300 
    301302
    302303#  Define center of the tree position inside the local LAD domain
     
    468469
    469470
    470 def process_patch(dz,patch_height,patch_lai,alpha,beta):
     471def process_patch(dz,patch_height,max_height_lad,patch_lai,alpha,beta):
    471472   
    472473#  Define fill values
     
    481482   pch_index = np.where( (pch_index[:,:] == -1) ,0,pch_index[:,:])
    482483
    483    max_canopy_height = max(patch_height.flatten())
     484   max_canopy_height = max(max(patch_height.flatten()),max_height_lad)
    484485
    485486   z = np.arange(0,math.floor(max_canopy_height/dz)*dz+2*dz,dz)
     
    500501          int_bpdf = 0.0
    501502          if ( (patch_height[j,i] != fillvalues["tree_data"]) & (patch_height[j,i] >= (0.5*dz)) ):
    502              for k in range(0,pch_index[j,i]):
     503             for k in range(1,pch_index[j,i]):
    503504                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] ) )
    504505
    505              for k in range(0,pch_index[j,i]):
     506             for k in range(1,pch_index[j,i]):
    506507                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]
    507508
     
    538539# dbh: default trunk diameter at breast height (1.4 m) (m)
    539540#
    540 #
    541541class tree:
    542542    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

    r3629 r3668  
    2525# -----------------
    2626# $Id$
     27# Some improvements and new routines
     28#
     29# 3629 2018-12-13 12:18:54Z maronga
    2730# Some new routines
    2831#
     
    3942
    4043from netCDF4 import Dataset
     44
     45def nc_read_from_file_1d_all(filename, varname):
     46
     47
     48   import numpy as np
     49   import sys
     50   
     51   try:
     52      f = open(filename)
     53      f.close()
     54#      print("Load: " + filename + ".")
     55   except FileNotFoundError:
     56      print("Error: " + filename + ". No such file. Aborting...")
     57      sys.exit(1)
     58   
     59   nc_file = Dataset(filename, "r+", format="NETCDF4")
     60   tmp_array = np.array(nc_file.variables[varname][:] , dtype=type(nc_file.variables[varname]))
     61
     62   return tmp_array
    4163
    4264def nc_read_from_file_2d(filename, varname, x0, x1, y0, y1):
     
    174196   return 0
    175197
     198def nc_write_global_attribute(filename,attribute,value):
     199
     200   import datetime
     201   import os
     202
     203   print("Writing attribute " + attribute + " to file...")
     204
     205   f =  Dataset(filename, "a", format="NETCDF4")
     206   
     207   f.setncattr(attribute,value)
     208
     209   f.close()
     210
    176211
    177212def nc_write_dimension(filename,varname,array,datatype):
  • palm/trunk/SCRIPTS/palm_csd_files/palm_csd_tools.py

    r3567 r3668  
    2525# -----------------
    2626# $Id$
     27# Minor changes
     28#
     29# 3567 2018-11-27 13:59:21Z maronga
    2730# Initial revisions
    28 #
    29 #
    30 #
    3131#
    3232#
     
    7777
    7878def interpolate_2d(array,x1,y1,x2,y2):
    79    
    80       tmp_int2d = interp2d(y1,x1,array,kind='linear')
    81       array_ip = tmp_int2d(y2.astype(float), x2.astype(float))
     79   tmp_int2d = interp2d(x1,y1,array,kind='linear')
     80   array_ip = tmp_int2d(x2.astype(float), y2.astype(float))
    8281
    83       return array_ip
     82   return array_ip
    8483   
    8584   
     
    136135   return result
    137136
     137def check_consistency_3(array1,array2,array3,fill1,fill2,fill3):
     138
     139   tmp_array = np.where(array1 != fill1,1,0) + np.where(array2 != fill2,1,0) + np.where(array3 != fill3,-1,0)
     140   
     141   test = np.any(tmp_array.flatten() != 0)
     142   if test:
     143      print("soil_type array is not consistent!")
     144      print("max: " + str(max(tmp_array.flatten())) + ", min: " + str(min(tmp_array.flatten())))
     145   else:
     146      print("soil_type array is consistent!")     
     147   return tmp_array, test
     148
     149
     150
    138151
    139152def check_consistency_4(array1,array2,array3,array4,fill1,fill2,fill3,fill4):
  • palm/trunk/SOURCE/check_parameters.f90

    r3655 r3668  
    2525! -----------------
    2626! $Id$
     27! most_method removed
     28!
     29! 3655 2019-01-07 16:51:22Z knoop
    2730! Formatting
    2831!
     
    37903793    ENDIF
    37913794
    3792 !
    3793 !-- Check for valid setting of most_method
    3794     IF ( TRIM( most_method ) /= 'circular'  .AND.                              &
    3795          TRIM( most_method ) /= 'newton'    .AND.                              &
    3796          TRIM( most_method ) /= 'lookup' )  THEN
    3797        message_string = 'most_method = "' // TRIM( most_method ) //            &
    3798                         '" is unknown'
    3799        CALL message( 'check_parameters', 'PA0416', 1, 2, 0, 6, 0 )
    3800     ENDIF
    38013795
    38023796!
  • palm/trunk/SOURCE/land_surface_model_mod.f90

    r3655 r3668  
    2525! -----------------
    2626! $Id$
     27! Removed most_method
     28!
     29! 3655 2019-01-07 16:51:22Z knoop
    2730! nopointer option removed
    2831!
     
    13241327
    13251328    USE control_parameters,                                                    &
    1326         ONLY:  bc_pt_b, bc_q_b, constant_flux_layer, message_string,           &
    1327                most_method
     1329        ONLY:  bc_pt_b, bc_q_b, constant_flux_layer, message_string
    13281330                     
    13291331   
     
    14461448   
    14471449    IF ( TRIM( surface_type ) == 'water' )  THEN
    1448 
    1449        IF ( TRIM( most_method ) == 'lookup' )  THEN   
    1450           WRITE( message_string, * ) 'surface_type = ', surface_type,          &
    1451                                      ' is not allowed in combination with ',   &
    1452                                      'most_method = ', most_method
    1453           CALL message( 'lsm_check_parameters', 'PA0414', 1, 2, 0, 6, 0 )
    1454        ENDIF
    14551450
    14561451       IF ( water_type == 0 )  THEN 
     
    15311526
    15321527    IF ( TRIM( surface_type ) == 'netcdf' )  THEN
    1533        IF ( ANY( water_type_f%var /= water_type_f%fill )  .AND.                &
    1534             TRIM( most_method ) == 'lookup' )  THEN   
    1535           WRITE( message_string, * ) 'water-surfaces are not allowed in ' //   &
    1536                                      'combination with most_method = ',        &
    1537                                      TRIM( most_method )
    1538           CALL message( 'lsm_check_parameters', 'PA0999', 2, 2, 0, 6, 0 )
    1539        ENDIF
    15401528!
    15411529!--    MS: Some problme here, after calling message everythings stucks at
  • palm/trunk/SOURCE/modules.f90

    r3648 r3668  
    2525! -----------------
    2626! $Id$
     27! Removed most_method
     28!
     29! 3648 2019-01-02 16:35:46Z suehring
    2730! -surface_data_output +surface_output
    2831!
     
    10901093    CHARACTER (LEN=8)    ::  coupling_char = ''                           !< appended to filenames in coupled or nested runs ('_O': ocean PE,
    10911094                                                                          !< '_NV': vertically nested atmosphere PE, '_N##': PE of nested domain ##
    1092     CHARACTER (LEN=8)    ::  most_method = 'newton'                       !< namelist parameter
    10931095    CHARACTER (LEN=10)   ::  run_date = ' '                               !< date of simulation run
    10941096    CHARACTER (LEN=8)    ::  run_time = ' '                               !< time of simulation run
  • palm/trunk/SOURCE/parin.f90

    r3649 r3668  
    2525! -----------------
    2626! $Id$
     27! Removed most_method
     28!
     29! 3649 2019-01-02 16:52:21Z suehring
    2730! Delete debug-print statements
    2831!
     
    573576             loop_optimization, lsf_exception, masking_method, mg_cycles,      &
    574577             mg_switch_to_pe0_level, mixing_length_1d, momentum_advec,         &
    575              most_method,                                                      &
    576578             netcdf_precision, neutral, ngsrb,                                 &
    577579             nsor, nsor_ini, nudging, nx, ny, nz, ocean_mode, omega,           &
     
    644646             loop_optimization, lsf_exception, masking_method, mg_cycles,      &
    645647             mg_switch_to_pe0_level, mixing_length_1d, momentum_advec,         &
    646              most_method,                                                      &
    647648             netcdf_precision, neutral, ngsrb,                                 &
    648649             nsor, nsor_ini, nudging, nx, ny, nz, ocean_mode, omega,           &
  • palm/trunk/SOURCE/read_restart_data_mod.f90

    r3655 r3668  
    2525! -----------------
    2626! $Id$
     27! Removed most_method and increased binary version
     28!
     29! 3655 2019-01-07 16:51:22Z knoop
    2730! Implementation of the PALM module interface
    2831!
     
    203206       READ ( 13 )  version_on_file
    204207
    205        binary_version_global = '4.7'
     208       binary_version_global = '4.8'
    206209       IF ( TRIM( version_on_file ) /= TRIM( binary_version_global ) )  THEN
    207210          WRITE( message_string, * ) 'version mismatch concerning ',           &
     
    488491             CASE ( 'momentum_advec' )
    489492                READ ( 13 )  momentum_advec
    490              CASE ( 'most_method' )
    491                 READ ( 13 )  most_method
    492493             CASE ( 'netcdf_precision' )
    493494                READ ( 13 )  netcdf_precision
  • palm/trunk/SOURCE/surface_layer_fluxes_mod.f90

    r3655 r3668  
    2626! -----------------
    2727! $Id$
     28! Removed methods "circular" and "lookup"
     29!
     30! 3655 2019-01-07 16:51:22Z knoop
    2831! OpenACC port for SPEC
    2932!
     
    228231! ------------
    229232!> Diagnostic computation of vertical fluxes in the constant flux layer from the
    230 !> values of the variables at grid point k=1. Three different methods are
    231 !> available:
    232 !> 1) the "old" version (most_method = 'circular') which is fast, but inaccurate
    233 !> 2) a Newton iteration method (most_method = 'newton'), which is accurate, but
    234 !>    slower
    235 !> 3) a method using a lookup table which is fast and accurate. Note, however,
    236 !>    that this method cannot be used in case of roughness heterogeneity
     233!> values of the variables at grid point k=1 based on Newton iteration
    237234!>
    238235!> @todo (re)move large_scale_forcing actions
     
    261258               intermediate_timestep_count, intermediate_timestep_count_max,   &
    262259               land_surface, large_scale_forcing, lsf_surf, message_string,    &
    263                most_method, neutral, passive_scalar, pt_surface, q_surface,    &
     260               neutral, passive_scalar, pt_surface, q_surface,                 &
    264261               run_coupled, surface_pressure, simulated_time, terminate_run,   &
    265262               time_since_reference_point, urban_surface,                      &
     
    400397!--    First, calculate the new Obukhov length, then new friction velocity,
    401398!--    followed by the new scaling parameters (th*, q*, etc.), and the new
    402 !--    surface fluxes if required. The old routine ("circular") requires a
    403 !--    different order of calls as the scaling parameters from the previous time
    404 !--    steps are used to calculate the Obukhov length
    405 
    406 !
    407 !--    Depending on setting of most_method use the "old" routine
    408 !--    Note, each routine is called for different surface types.
     399!--    surface fluxes if required. Note, each routine is called for different surface types.
    409400!--    First call for default-type horizontal surfaces, for natural- and
    410401!--    urban-type surfaces. Note, at this place only upward-facing horizontal
    411 !--    surfaces are treted.
    412        IF ( most_method == 'circular' )  THEN
    413 !
    414 !--       Default-type upward-facing horizontal surfaces
    415           IF ( surf_def_h(0)%ns >= 1  )  THEN
    416              surf => surf_def_h(0)
    417              CALL calc_scaling_parameters
    418              CALL calc_uvw_abs
    419              IF ( .NOT. neutral )  CALL calc_ol
    420              CALL calc_us
    421              CALL calc_surface_fluxes
    422              
    423              IF ( do_output_at_2m )  THEN
    424                 CALL calc_pt_near_surface ( '2m' )
    425              ENDIF
    426           ENDIF
    427 !
    428 !--       Natural-type horizontal surfaces
    429           IF ( surf_lsm_h%ns >= 1  )  THEN
    430              surf => surf_lsm_h
    431              CALL calc_scaling_parameters
    432              CALL calc_uvw_abs
    433              IF ( .NOT. neutral )  CALL calc_ol
    434              CALL calc_us
    435              CALL calc_surface_fluxes
    436              IF ( do_output_at_2m )  THEN
    437                 CALL calc_pt_near_surface ( '2m' )
    438              ENDIF
    439           ENDIF
    440 !
    441 !--       Urban-type horizontal surfaces
    442           IF ( surf_usm_h%ns >= 1  )  THEN
    443              surf => surf_usm_h
    444              CALL calc_scaling_parameters
    445              CALL calc_uvw_abs
    446              IF ( .NOT. neutral )  CALL calc_ol
    447              CALL calc_us
    448              CALL calc_surface_fluxes
    449              IF ( do_output_at_2m )  THEN
    450                 CALL calc_pt_near_surface ( '2m' )
    451              ENDIF
    452              IF ( indoor_model )  THEN
    453                 CALL calc_pt_near_surface ( '10cm' )
    454              ENDIF
    455           ENDIF
    456 !
    457 !--    Use either Newton iteration or a lookup table for the bulk Richardson
    458 !--    number to calculate the Obukhov length
    459        ELSEIF ( most_method == 'newton'  .OR.  most_method == 'lookup' )  THEN
    460 !
    461 !--       Default-type upward-facing horizontal surfaces
    462           IF ( surf_def_h(0)%ns >= 1  )  THEN
    463              surf => surf_def_h(0)
    464              CALL calc_uvw_abs
    465              IF ( .NOT. neutral )  CALL calc_ol
    466              CALL calc_us
    467              CALL calc_scaling_parameters
    468              CALL calc_surface_fluxes
    469              IF ( do_output_at_2m )  THEN
    470                 CALL calc_pt_near_surface ( '2m' )
    471              ENDIF
    472           ENDIF
    473 !
    474 !--       Natural-type horizontal surfaces
    475           IF ( surf_lsm_h%ns >= 1  )  THEN
    476              surf => surf_lsm_h
    477              CALL calc_uvw_abs
    478              IF ( .NOT. neutral )  CALL calc_ol
    479              CALL calc_us
    480              CALL calc_scaling_parameters
    481              CALL calc_surface_fluxes
    482              IF ( do_output_at_2m )  THEN
    483                 CALL calc_pt_near_surface ( '2m' )
    484              ENDIF
    485           ENDIF
    486 !
    487 !--       Urban-type horizontal surfaces
    488           IF ( surf_usm_h%ns >= 1  )  THEN
    489              surf => surf_usm_h
    490              CALL calc_uvw_abs
    491              IF ( .NOT. neutral )  CALL calc_ol
    492              CALL calc_us
    493              CALL calc_scaling_parameters
    494              CALL calc_surface_fluxes
    495              IF ( do_output_at_2m )  THEN
    496                 CALL calc_pt_near_surface ( '2m' )
    497              ENDIF
    498              IF ( indoor_model )  THEN
    499                 CALL calc_pt_near_surface ( '10cm' )
    500              ENDIF
    501           ENDIF
    502 
     402!--    surfaces are treated.
     403
     404!
     405!--    Default-type upward-facing horizontal surfaces
     406       IF ( surf_def_h(0)%ns >= 1  )  THEN
     407          surf => surf_def_h(0)
     408          CALL calc_uvw_abs
     409          IF ( .NOT. neutral )  CALL calc_ol
     410          CALL calc_us
     411          CALL calc_scaling_parameters
     412          CALL calc_surface_fluxes
     413          IF ( do_output_at_2m )  THEN
     414             CALL calc_pt_near_surface ( '2m' )
     415          ENDIF
    503416       ENDIF
     417!
     418!--    Natural-type horizontal surfaces
     419       IF ( surf_lsm_h%ns >= 1  )  THEN
     420          surf => surf_lsm_h
     421          CALL calc_uvw_abs
     422          IF ( .NOT. neutral )  CALL calc_ol
     423          CALL calc_us
     424          CALL calc_scaling_parameters
     425          CALL calc_surface_fluxes
     426          IF ( do_output_at_2m )  THEN
     427             CALL calc_pt_near_surface ( '2m' )
     428          ENDIF
     429       ENDIF
     430!
     431!--    Urban-type horizontal surfaces
     432       IF ( surf_usm_h%ns >= 1  )  THEN
     433          surf => surf_usm_h
     434          CALL calc_uvw_abs
     435          IF ( .NOT. neutral )  CALL calc_ol
     436          CALL calc_us
     437          CALL calc_scaling_parameters
     438          CALL calc_surface_fluxes
     439          IF ( do_output_at_2m )  THEN
     440             CALL calc_pt_near_surface ( '2m' )
     441          ENDIF
     442          IF ( indoor_model )  THEN
     443             CALL calc_pt_near_surface ( '10cm' )
     444          ENDIF
     445       ENDIF
     446
    504447!
    505448!--    Treat downward-facing horizontal surfaces. Note, so far, these are
     
    549492!--    flux in land-surface model in case of aero_resist_kray is not true.
    550493       IF ( .NOT. aero_resist_kray )  THEN
    551           IF ( most_method == 'circular' )  THEN
    552              DO  l = 0, 1
    553                 IF ( surf_lsm_v(l)%ns >= 1  )  THEN
    554                    surf => surf_lsm_v(l)
    555 !
    556 !--                Compute scaling parameters
    557                    CALL calc_scaling_parameters
    558 !
    559 !--                Compute surface-parallel velocity
    560                    CALL calc_uvw_abs_v_ugrid
    561 !
    562 !--                Compute Obukhov length
    563                    IF ( .NOT. neutral )  CALL calc_ol
    564 !
    565 !--                Compute respective friction velocity on staggered grid
    566                    CALL calc_us
    567 !
    568 !--                Compute respective surface fluxes for momentum and TKE
    569                    CALL calc_surface_fluxes
    570                 ENDIF
    571              ENDDO
    572           ELSE
    573              DO  l = 0, 1
    574                 IF ( surf_lsm_v(l)%ns >= 1  )  THEN
    575                    surf => surf_lsm_v(l)
    576 !
    577 !--                Compute surface-parallel velocity
    578                    CALL calc_uvw_abs_v_ugrid
    579 !
    580 !--                Compute Obukhov length
    581                    IF ( .NOT. neutral )  CALL calc_ol
    582 !
    583 !--                Compute respective friction velocity on staggered grid
    584                    CALL calc_us
    585 !
    586 !--                Compute scaling parameters
    587                    CALL calc_scaling_parameters
    588 !
    589 !--                Compute respective surface fluxes for momentum and TKE
    590                    CALL calc_surface_fluxes
    591                 ENDIF
    592              ENDDO
    593           ENDIF
     494          DO  l = 0, 1
     495             IF ( surf_lsm_v(l)%ns >= 1  )  THEN
     496                surf => surf_lsm_v(l)
     497!
     498!--             Compute surface-parallel velocity
     499                CALL calc_uvw_abs_v_ugrid
     500!
     501!--             Compute Obukhov length
     502                IF ( .NOT. neutral )  CALL calc_ol
     503!
     504!--             Compute respective friction velocity on staggered grid
     505                CALL calc_us
     506!
     507!--             Compute scaling parameters
     508                CALL calc_scaling_parameters
     509!
     510!--             Compute respective surface fluxes for momentum and TKE
     511                CALL calc_surface_fluxes
     512             ENDIF
     513          ENDDO
    594514!
    595515!--    No ts is required, so scaling parameters and Obukhov length do not need
     
    652572!--    flux in land-surface model in case of aero_resist_kray is not true.
    653573       IF ( .NOT. aero_resist_kray )  THEN
    654           IF ( most_method == 'circular' )  THEN
    655              DO  l = 2, 3
    656                 IF ( surf_lsm_v(l)%ns >= 1  )  THEN
    657                    surf => surf_lsm_v(l)
    658 !
    659 !--                Compute scaling parameters
    660                    CALL calc_scaling_parameters
    661 !
    662 !--                Compute surface-parallel velocity
    663                    CALL calc_uvw_abs_v_vgrid
    664 !
    665 !--                Compute Obukhov length
    666                    IF ( .NOT. neutral )  CALL calc_ol
    667 !
    668 !--                Compute respective friction velocity on staggered grid
    669                    CALL calc_us
    670 !
    671 !--                Compute respective surface fluxes for momentum and TKE
    672                    CALL calc_surface_fluxes
    673                 ENDIF
    674              ENDDO
    675           ELSE
    676              DO  l = 2, 3
    677                 IF ( surf_lsm_v(l)%ns >= 1  )  THEN
    678                    surf => surf_lsm_v(l)
    679 !
    680 !--                Compute surface-parallel velocity
    681                    CALL calc_uvw_abs_v_vgrid
    682 !
    683 !--                Compute Obukhov length
    684                    IF ( .NOT. neutral )  CALL calc_ol
    685 !
    686 !--                Compute respective friction velocity on staggered grid
    687                    CALL calc_us
    688 !
    689 !--                Compute scaling parameters
    690                    CALL calc_scaling_parameters
    691 !
    692 !--                Compute respective surface fluxes for momentum and TKE
    693                    CALL calc_surface_fluxes
    694                 ENDIF
    695              ENDDO
    696           ENDIF
     574          DO  l = 2, 3
     575             IF ( surf_lsm_v(l)%ns >= 1  )  THEN
     576                surf => surf_lsm_v(l)
     577!
     578!--             Compute surface-parallel velocity
     579                CALL calc_uvw_abs_v_vgrid
     580!
     581!--             Compute Obukhov length
     582                IF ( .NOT. neutral )  CALL calc_ol
     583!
     584!--             Compute respective friction velocity on staggered grid
     585                CALL calc_us
     586!
     587!--             Compute scaling parameters
     588                CALL calc_scaling_parameters
     589!
     590!--             Compute respective surface fluxes for momentum and TKE
     591                CALL calc_surface_fluxes
     592             ENDIF
     593          ENDDO
    697594       ELSE
    698595          DO  l = 2, 3
     
    841738          IF ( surf_lsm_h%ns    >= 1 )  surf_lsm_h%ol    = 1.0E10_wp
    842739          IF ( surf_usm_h%ns    >= 1 )  surf_usm_h%ol    = 1.0E10_wp
    843        ENDIF
    844 
    845        IF ( most_method == 'lookup' )  THEN
    846 
    847 !
    848 !--       Check for roughness heterogeneity. In that case terminate run and
    849 !--       inform user. Check for both, natural and non-natural walls.
    850           IF ( surf_def_h(0)%ns >= 1 )  THEN
    851              IF ( MINVAL( surf_def_h(0)%z0h ) /= MAXVAL( surf_def_h(0)%z0h )  .OR. &
    852                   MINVAL( surf_def_h(0)%z0  ) /= MAXVAL( surf_def_h(0)%z0  ) )  THEN
    853                 terminate_run_l = .TRUE.
    854              ENDIF
    855           ENDIF
    856           IF ( surf_lsm_h%ns >= 1 )  THEN
    857              IF ( MINVAL( surf_lsm_h%z0h ) /= MAXVAL( surf_lsm_h%z0h )  .OR.       &
    858                   MINVAL( surf_lsm_h%z0  ) /= MAXVAL( surf_lsm_h%z0  ) )  THEN
    859                 terminate_run_l = .TRUE.
    860              ENDIF
    861           ENDIF
    862           IF ( surf_usm_h%ns >= 1 )  THEN
    863              IF ( MINVAL( surf_usm_h%z0h ) /= MAXVAL( surf_usm_h%z0h )  .OR.       &
    864                   MINVAL( surf_usm_h%z0  ) /= MAXVAL( surf_usm_h%z0  ) )  THEN
    865                 terminate_run_l = .TRUE.
    866              ENDIF
    867           ENDIF
    868 !
    869 !--       Check roughness homogeneity between differt surface types.
    870           IF ( surf_lsm_h%ns >= 1  .AND.  surf_def_h(0)%ns >= 1 )  THEN
    871              IF ( MINVAL( surf_lsm_h%z0h ) /= MAXVAL( surf_def_h(0)%z0h )  .OR.    &
    872                   MINVAL( surf_lsm_h%z0  ) /= MAXVAL( surf_def_h(0)%z0  ) )  THEN
    873                 terminate_run_l = .TRUE.
    874              ENDIF
    875           ENDIF
    876           IF ( surf_usm_h%ns >= 1  .AND.  surf_def_h(0)%ns >= 1 )  THEN
    877              IF ( MINVAL( surf_usm_h%z0h ) /= MAXVAL( surf_def_h(0)%z0h )  .OR.    &
    878                   MINVAL( surf_usm_h%z0  ) /= MAXVAL( surf_def_h(0)%z0  ) )  THEN
    879                 terminate_run_l = .TRUE.
    880              ENDIF
    881           ENDIF
    882           IF ( surf_usm_h%ns >= 1  .AND.  surf_lsm_h%ns >= 1 )  THEN
    883              IF ( MINVAL( surf_usm_h%z0h ) /= MAXVAL( surf_lsm_h%z0h )  .OR.       &
    884                   MINVAL( surf_usm_h%z0  ) /= MAXVAL( surf_lsm_h%z0  ) )  THEN
    885                 terminate_run_l = .TRUE.
    886              ENDIF
    887           ENDIF
    888 
    889 
    890 #if defined( __parallel )
    891 !
    892 !--       Make a logical OR for all processes. Force termiation of model if result
    893 !--       is TRUE
    894           IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    895           CALL MPI_ALLREDUCE( terminate_run_l, terminate_run, 1, MPI_LOGICAL,  &
    896                               MPI_LOR, comm2d, ierr )
    897 #else
    898           terminate_run = terminate_run_l
    899 #endif
    900 
    901           IF ( terminate_run )  THEN
    902              message_string = 'most_method = "lookup" cannot be used in ' //   &
    903                               'combination with a prescribed roughness '  //   &
    904                               'heterogeneity'
    905              CALL message( 'surface_layer_fluxes', 'PA0116', 1, 2, 0, 6, 0 )
    906           ENDIF
    907 
    908           ALLOCATE(  zeta_tmp(0:num_steps-1) )
    909 
    910 !
    911 !--       Use the lowest possible value for z_mo
    912           k    = nzb
    913           z_mo = zu(k+1) - zw(k)
    914 
    915 !
    916 !--       Calculate z/L range from zeta_stretch to zeta_max using 90% of the
    917 !--       available steps (num_steps). The calculation is done with negative
    918 !--       values of zeta in order to simplify the stretching in the free
    919 !--       convection limit for the remaining 10% of steps.
    920           zeta_tmp(0) = - zeta_max
    921           num_steps_n = ( num_steps * 9 / 10 ) - 1
    922           zeta_step   = (zeta_max - zeta_stretch) / REAL(num_steps_n)
    923 
    924           DO li = 1, num_steps_n
    925              zeta_tmp(li) = zeta_tmp(li-1) + zeta_step
    926           ENDDO
    927 
    928 !
    929 !--       Calculate stretching factor for the free convection range
    930           DO  WHILE ( ABS( (regr-regr_old) / regr_old ) > 1.0E-10_wp )
    931              regr_old = regr
    932              regr = ( 1.0_wp - ( -zeta_min / zeta_step ) * ( 1.0_wp - regr )   &
    933                     )**( 10.0_wp / REAL(num_steps) )
    934           ENDDO
    935 
    936 !
    937 !--       Calculate z/L range from zeta_min to zeta_stretch
    938           DO li = num_steps_n+1, num_steps-1
    939              zeta_tmp(li) = zeta_tmp(li-1) + zeta_step
    940              zeta_step = zeta_step * regr
    941           ENDDO
    942 
    943 !
    944 !--       Save roughness lengths to temporary variables
    945           IF ( surf_def_h(0)%ns >= 1  )  THEN
    946              z0h_min = surf_def_h(0)%z0h(1)
    947              z0_min  = surf_def_h(0)%z0(1)
    948           ELSEIF ( surf_lsm_h%ns >= 1 )  THEN
    949              z0h_min = surf_lsm_h%z0h(1)
    950              z0_min  = surf_lsm_h%z0(1)
    951           ELSEIF ( surf_usm_h%ns >= 1 )  THEN
    952              z0h_min = surf_usm_h%z0h(1)
    953              z0_min  = surf_usm_h%z0(1)
    954           ENDIF         
    955 !
    956 !--       Calculate lookup table for the Richardson number versus Obukhov length
    957 !--       The Richardson number (rib) is defined depending on the choice of
    958 !--       boundary conditions for temperature
    959           IF ( ibc_pt_b == 1 )  THEN
    960              DO li = 0, num_steps-1
    961                 ol_tab(li)  = - z_mo / zeta_tmp(num_steps-1-li)
    962                 rib_tab(li) = z_mo / ol_tab(li)  / ( LOG( z_mo / z0_min )       &
    963                                                 - psi_m( z_mo / ol_tab(li) )    &
    964                                                 + psi_m( z0_min / ol_tab(li) )  &
    965                                                   )**3
    966              ENDDO 
    967           ELSE
    968              DO li = 0, num_steps-1
    969                 ol_tab(li)  = - z_mo / zeta_tmp(num_steps-1-li)
    970                 rib_tab(li) = z_mo / ol_tab(li)  * ( LOG( z_mo / z0h_min )     &
    971                                               - psi_h( z_mo / ol_tab(li) )     &
    972                                               + psi_h( z0h_min / ol_tab(li) )  &
    973                                             )                                  &
    974                                           / ( LOG( z_mo / z0_min )             &
    975                                               - psi_m( z_mo / ol_tab(li) )     &
    976                                               + psi_m( z0_min / ol_tab(li) )   &
    977                                             )**2
    978              ENDDO
    979           ENDIF
    980 
    981 !
    982 !--       Determine minimum values of rib in the lookup table. Set upper limit
    983 !--       to critical Richardson number (0.25)
    984           rib_min  = MINVAL(rib_tab)
    985           rib_max  = 0.25 !MAXVAL(rib_tab)
    986 
    987           DEALLOCATE( zeta_tmp )
    988740       ENDIF
    989741
     
    12561008                       ol_u      !< Upper bound of L for Newton iteration
    12571009
    1258        IF ( TRIM( most_method ) /= 'circular' )  THEN
    1259 !
    1260 !--       Evaluate bulk Richardson number (calculation depends on
    1261 !--       definition based on setting of boundary conditions
    1262           IF ( ibc_pt_b /= 1 )  THEN
    1263              IF ( humidity )  THEN
    1264                 !$OMP PARALLEL DO PRIVATE( z_mo )
    1265                 DO  m = 1, surf%ns
    1266 
    1267                    z_mo = surf%z_mo(m)
    1268 
    1269                    surf%rib(m) = g * z_mo                                      &
    1270                                  * ( surf%vpt1(m) - surf%vpt_surface(m) )      &
    1271                                  / ( surf%uvw_abs(m)**2 * surf%vpt1(m)         &
    1272                                  + 1.0E-20_wp )
    1273                 ENDDO
    1274              ELSE
    1275                 !$OMP PARALLEL DO PRIVATE( z_mo )
    1276                 DO  m = 1, surf%ns
    1277 
    1278                    z_mo = surf%z_mo(m)
    1279 
    1280                    surf%rib(m) = g * z_mo                                      &
    1281                                  * ( surf%pt1(m) - surf%pt_surface(m) )        &
    1282                                  / ( surf%uvw_abs(m)**2 * surf%pt1(m)          &
    1283                                  + 1.0E-20_wp )
    1284                 ENDDO
    1285              ENDIF
     1010!
     1011!--    Evaluate bulk Richardson number (calculation depends on
     1012!--    definition based on setting of boundary conditions
     1013       IF ( ibc_pt_b /= 1 )  THEN
     1014          IF ( humidity )  THEN
     1015             !$OMP PARALLEL DO PRIVATE( z_mo )
     1016             DO  m = 1, surf%ns
     1017
     1018                z_mo = surf%z_mo(m)
     1019
     1020                surf%rib(m) = g * z_mo                                         &
     1021                              * ( surf%vpt1(m) - surf%vpt_surface(m) )         &
     1022                              / ( surf%uvw_abs(m)**2 * surf%vpt1(m)            &
     1023                              + 1.0E-20_wp )
     1024             ENDDO
    12861025          ELSE
    1287              IF ( humidity )  THEN
    1288                 !$OMP PARALLEL DO PRIVATE( k, z_mo )
    1289                 DO  m = 1, surf%ns
    1290 
    1291                    k   = surf%k(m)
    1292 
    1293                    z_mo = surf%z_mo(m)
    1294 
    1295                    surf%rib(m) = - g * z_mo * ( ( 1.0_wp + 0.61_wp             &
     1026             !$OMP PARALLEL DO PRIVATE( z_mo )
     1027             DO  m = 1, surf%ns
     1028
     1029                z_mo = surf%z_mo(m)
     1030
     1031                surf%rib(m) = g * z_mo                                         &
     1032                              * ( surf%pt1(m) - surf%pt_surface(m) )           &
     1033                              / ( surf%uvw_abs(m)**2 * surf%pt1(m) + 1.0E-20_wp )
     1034             ENDDO
     1035          ENDIF
     1036       ELSE
     1037          IF ( humidity )  THEN
     1038             !$OMP PARALLEL DO PRIVATE( k, z_mo )
     1039             DO  m = 1, surf%ns
     1040
     1041                k   = surf%k(m)
     1042
     1043                z_mo = surf%z_mo(m)
     1044
     1045                surf%rib(m) = - g * z_mo * ( ( 1.0_wp + 0.61_wp                &
    12961046                           * surf%qv1(m) ) * surf%shf(m) + 0.61_wp             &
    12971047                           * surf%pt1(m) * surf%qsws(m) ) *                    &
     
    12991049                         ( surf%uvw_abs(m)**3 * surf%vpt1(m) * kappa**2        &
    13001050                           + 1.0E-20_wp )
    1301                 ENDDO
     1051             ENDDO
     1052          ELSE
     1053             !$OMP PARALLEL DO PRIVATE( k, z_mo )
     1054             !$ACC PARALLEL LOOP PRIVATE(k, z_mo) &
     1055             !$ACC PRESENT(surf, drho_air_zw)
     1056             DO  m = 1, surf%ns
     1057
     1058                k   = surf%k(m)
     1059
     1060                z_mo = surf%z_mo(m)
     1061
     1062                surf%rib(m) = - g * z_mo * surf%shf(m) *                    &
     1063                                     drho_air_zw(k-1)            /          &
     1064                        ( surf%uvw_abs(m)**3 * surf%pt1(m) * kappa**2       &
     1065                        + 1.0E-20_wp )
     1066             ENDDO
     1067          ENDIF
     1068       ENDIF
     1069
     1070!
     1071!--    Calculate the Obukhov length using Newton iteration
     1072       !$ACC PARALLEL LOOP PRIVATE(i, j, z_mo) &
     1073       !$ACC PRIVATE(ol_old, ol_m, ol_l, ol_u, f, f_d_ol) &
     1074       !$ACC PRESENT(surf)
     1075       DO  m = 1, surf%ns
     1076
     1077          i   = surf%i(m)           
     1078          j   = surf%j(m)
     1079
     1080          z_mo = surf%z_mo(m)
     1081
     1082!
     1083!--       Store current value in case the Newton iteration fails
     1084          ol_old = surf%ol(m)
     1085
     1086!
     1087!--       Ensure that the bulk Richardson number and the Obukhov
     1088!--       length have the same sign
     1089          IF ( surf%rib(m) * surf%ol(m) < 0.0_wp  .OR.                      &
     1090               ABS( surf%ol(m) ) == ol_max )  THEN
     1091             IF ( surf%rib(m) > 1.0_wp ) surf%ol(m) =  0.01_wp
     1092             IF ( surf%rib(m) < 0.0_wp ) surf%ol(m) = -0.01_wp
     1093          ENDIF
     1094!
     1095!--       Iteration to find Obukhov length
     1096          iter = 0
     1097          DO
     1098             iter = iter + 1
     1099!
     1100!--          In case of divergence, use the value of the previous time step
     1101             IF ( iter > 1000 )  THEN
     1102                surf%ol(m) = ol_old
     1103                EXIT
     1104             ENDIF
     1105
     1106             ol_m = surf%ol(m)
     1107             ol_l = ol_m - 0.001_wp * ol_m
     1108             ol_u = ol_m + 0.001_wp * ol_m
     1109
     1110
     1111             IF ( ibc_pt_b /= 1 )  THEN
     1112!
     1113!--             Calculate f = Ri - [...]/[...]^2 = 0
     1114                f = surf%rib(m) - ( z_mo / ol_m ) * (                          &
     1115                                              LOG( z_mo / surf%z0h(m) )        &
     1116                                              - psi_h( z_mo / ol_m )           &
     1117                                              + psi_h( surf%z0h(m) /           &
     1118                                                       ol_m )                  & 
     1119                                                     )                         &
     1120                                           / ( LOG( z_mo / surf%z0(m) )        &
     1121                                              - psi_m( z_mo / ol_m )           &
     1122                                              + psi_m( surf%z0(m) /  ol_m )    &
     1123                                                 )**2
     1124
     1125!
     1126!--              Calculate df/dL
     1127                 f_d_ol = ( - ( z_mo / ol_u ) * ( LOG( z_mo /                  &
     1128                                                          surf%z0h(m) )        &
     1129                                         - psi_h( z_mo / ol_u )                &
     1130                                         + psi_h( surf%z0h(m) / ol_u )         &
     1131                                           )                                   &
     1132                                         / ( LOG( z_mo / surf%z0(m) )          &
     1133                                         - psi_m( z_mo / ol_u )                &
     1134                                         + psi_m( surf%z0(m) / ol_u )          &
     1135                                           )**2                                &
     1136                        + ( z_mo / ol_l ) * ( LOG( z_mo / surf%z0h(m) )        &
     1137                                         - psi_h( z_mo / ol_l )                &
     1138                                         + psi_h( surf%z0h(m) / ol_l )         &
     1139                                           )                                   &
     1140                                         / ( LOG( z_mo / surf%z0(m) )          &
     1141                                         - psi_m( z_mo / ol_l )                &
     1142                                         + psi_m( surf%z0(m) / ol_l )          &
     1143                                           )**2                                &
     1144                          ) / ( ol_u - ol_l )
    13021145             ELSE
    1303                 !$OMP PARALLEL DO PRIVATE( k, z_mo )
    1304                 !$ACC PARALLEL LOOP PRIVATE(k, z_mo) &
    1305                 !$ACC PRESENT(surf, drho_air_zw)
    1306                 DO  m = 1, surf%ns
    1307 
    1308                    k   = surf%k(m)
    1309 
    1310                    z_mo = surf%z_mo(m)
    1311 
    1312                    surf%rib(m) = - g * z_mo * surf%shf(m) *                    &
    1313                                         drho_air_zw(k-1)            /          &
    1314                         ( surf%uvw_abs(m)**3 * surf%pt1(m) * kappa**2          &
    1315                            + 1.0E-20_wp )
    1316                 ENDDO
     1146!
     1147!--             Calculate f = Ri - 1 /[...]^3 = 0
     1148                f = surf%rib(m) - ( z_mo / ol_m ) /                            &
     1149                                             ( LOG( z_mo / surf%z0(m) )        &
     1150                                         - psi_m( z_mo / ol_m )                &
     1151                                         + psi_m( surf%z0(m) / ol_m )          &
     1152                                             )**3
     1153
     1154!
     1155!--             Calculate df/dL
     1156                f_d_ol = ( - ( z_mo / ol_u ) / ( LOG( z_mo / surf%z0(m) )      &
     1157                                         - psi_m( z_mo / ol_u )                &
     1158                                         + psi_m( surf%z0(m) / ol_u )          &
     1159                                                  )**3                         &
     1160                           + ( z_mo / ol_l ) / ( LOG( z_mo / surf%z0(m) )      &
     1161                                         - psi_m( z_mo / ol_l )                &
     1162                                         + psi_m( surf%z0(m) / ol_l )          &
     1163                                            )**3                               &
     1164                          ) / ( ol_u - ol_l )
    13171165             ENDIF
    1318           ENDIF
    1319 
    1320        ENDIF
    1321 
    1322 
    1323 !
    1324 !--    Calculate the Obukhov length either using a Newton iteration
    1325 !--    method, via a lookup table, or using the old circular way
    1326        IF ( TRIM( most_method ) == 'newton' )  THEN
    1327 
    1328           !$ACC PARALLEL LOOP PRIVATE(i, j, z_mo) &
    1329           !$ACC PRIVATE(ol_old, ol_m, ol_l, ol_u, f, f_d_ol) &
    1330           !$ACC PRESENT(surf)
    1331           DO  m = 1, surf%ns
    1332 
    1333              i   = surf%i(m)           
    1334              j   = surf%j(m)
    1335 
    1336              z_mo = surf%z_mo(m)
    1337 
    1338 !
    1339 !--          Store current value in case the Newton iteration fails
    1340              ol_old = surf%ol(m)
     1166!
     1167!--          Calculate new L
     1168             surf%ol(m) = ol_m - f / f_d_ol
    13411169
    13421170!
    13431171!--          Ensure that the bulk Richardson number and the Obukhov
    1344 !--          length have the same sign
    1345              IF ( surf%rib(m) * surf%ol(m) < 0.0_wp  .OR.                      &
    1346                   ABS( surf%ol(m) ) == ol_max )  THEN
    1347                 IF ( surf%rib(m) > 1.0_wp ) surf%ol(m) =  0.01_wp
    1348                 IF ( surf%rib(m) < 0.0_wp ) surf%ol(m) = -0.01_wp
     1172!--          length have the same sign and ensure convergence.
     1173             IF ( surf%ol(m) * ol_m < 0.0_wp )  surf%ol(m) = ol_m * 0.5_wp
     1174!
     1175!--          If unrealistic value occurs, set L to the maximum
     1176!--          value that is allowed
     1177             IF ( ABS( surf%ol(m) ) > ol_max )  THEN
     1178                surf%ol(m) = ol_max
     1179                EXIT
    13491180             ENDIF
    13501181!
    1351 !--          Iteration to find Obukhov length
    1352              iter = 0
    1353              DO
    1354                 iter = iter + 1
    1355 !
    1356 !--             In case of divergence, use the value of the previous time step
    1357                 IF ( iter > 1000 )  THEN
    1358                    surf%ol(m) = ol_old
    1359                    EXIT
    1360                 ENDIF
    1361 
    1362                 ol_m = surf%ol(m)
    1363                 ol_l = ol_m - 0.001_wp * ol_m
    1364                 ol_u = ol_m + 0.001_wp * ol_m
    1365 
    1366 
    1367                 IF ( ibc_pt_b /= 1 )  THEN
    1368 !
    1369 !--                Calculate f = Ri - [...]/[...]^2 = 0
    1370                    f = surf%rib(m) - ( z_mo / ol_m ) * (                       &
    1371                                                  LOG( z_mo / surf%z0h(m) )     &
    1372                                                  - psi_h( z_mo / ol_m )        &
    1373                                                  + psi_h( surf%z0h(m) /        &
    1374                                                           ol_m )               &
    1375                                                                )               &
    1376                                               / ( LOG( z_mo / surf%z0(m) )     &
    1377                                                  - psi_m( z_mo / ol_m )        &
    1378                                                  + psi_m( surf%z0(m) /         &
    1379                                                           ol_m )               &
    1380                                                  )**2
    1381 
    1382 !
    1383 !--                 Calculate df/dL
    1384                     f_d_ol = ( - ( z_mo / ol_u ) * ( LOG( z_mo /               &
    1385                                                              surf%z0h(m) )     &
    1386                                             - psi_h( z_mo / ol_u )             &
    1387                                             + psi_h( surf%z0h(m) / ol_u )      &
    1388                                               )                                &
    1389                                             / ( LOG( z_mo / surf%z0(m) )       &
    1390                                             - psi_m( z_mo / ol_u )             &
    1391                                             + psi_m( surf%z0(m) / ol_u )       &
    1392                                               )**2                             &
    1393                            + ( z_mo / ol_l ) * ( LOG( z_mo / surf%z0h(m) )     &
    1394                                             - psi_h( z_mo / ol_l )             &
    1395                                             + psi_h( surf%z0h(m) / ol_l )      &
    1396                                               )                                &
    1397                                             / ( LOG( z_mo / surf%z0(m) )       &
    1398                                             - psi_m( z_mo / ol_l )             &
    1399                                             + psi_m( surf%z0(m) / ol_l )       &
    1400                                               )**2                             &
    1401                              ) / ( ol_u - ol_l )
    1402                 ELSE
    1403 !
    1404 !--                Calculate f = Ri - 1 /[...]^3 = 0
    1405                    f = surf%rib(m) - ( z_mo / ol_m ) /                         &
    1406                                                 ( LOG( z_mo / surf%z0(m) )     &
    1407                                             - psi_m( z_mo / ol_m )             &
    1408                                             + psi_m( surf%z0(m) / ol_m )       &
    1409                                                 )**3
    1410 
    1411 !
    1412 !--                Calculate df/dL
    1413                    f_d_ol = ( - ( z_mo / ol_u ) / ( LOG( z_mo /                &
    1414                                                             surf%z0(m) )       &
    1415                                             - psi_m( z_mo / ol_u )             &
    1416                                             + psi_m( surf%z0(m) / ol_u )       &
    1417                                                      )**3                      &
    1418                            + ( z_mo / ol_l ) / ( LOG( z_mo / surf%z0(m) )      &
    1419                                             - psi_m( z_mo / ol_l )             &
    1420                                             + psi_m( surf%z0(m) / ol_l )       &
    1421                                                )**3                            &
    1422                                   ) / ( ol_u - ol_l )
    1423                 ENDIF
    1424 !
    1425 !--             Calculate new L
    1426                 surf%ol(m) = ol_m - f / f_d_ol
    1427 
    1428 !
    1429 !--             Ensure that the bulk Richardson number and the Obukhov
    1430 !--             length have the same sign and ensure convergence.
    1431                 IF ( surf%ol(m) * ol_m < 0.0_wp )  surf%ol(m) = ol_m * 0.5_wp
    1432 !
    1433 !--             If unrealistic value occurs, set L to the maximum
    1434 !--             value that is allowed
    1435                 IF ( ABS( surf%ol(m) ) > ol_max )  THEN
    1436                    surf%ol(m) = ol_max
    1437                    EXIT
    1438                 ENDIF
    1439 !
    1440 !--             Check for convergence
    1441                 IF ( ABS( ( surf%ol(m) - ol_m ) /                              &
    1442                      surf%ol(m) ) < 1.0E-4_wp )  THEN
    1443                    EXIT
    1444                 ELSE
    1445                    CYCLE
    1446                 ENDIF
    1447 
    1448              ENDDO
     1182!--          Check for convergence
     1183             IF ( ABS( ( surf%ol(m) - ol_m ) /  surf%ol(m) ) < 1.0E-4_wp )  THEN
     1184                EXIT
     1185             ELSE
     1186                CYCLE
     1187             ENDIF
     1188
    14491189          ENDDO
    1450 
    1451        ELSEIF ( TRIM( most_method ) == 'lookup' )  THEN
    1452 
    1453           !$OMP PARALLEL DO PRIVATE( i, j, z_mo, li ) FIRSTPRIVATE( li_bnd ) LASTPRIVATE( li_bnd )
    1454           DO  m = 1, surf%ns
    1455 
    1456              i   = surf%i(m)           
    1457              j   = surf%j(m)
    1458 !
    1459 !--          If the bulk Richardson number is outside the range of the lookup
    1460 !--          table, set it to the exceeding threshold value
    1461              IF ( surf%rib(m) < rib_min )  surf%rib(m) = rib_min
    1462              IF ( surf%rib(m) > rib_max )  surf%rib(m) = rib_max
    1463 
    1464 !
    1465 !--          Find the correct index bounds for linear interpolation. As the
    1466 !--          Richardson number will not differ very much from time step to
    1467 !--          time step , use the index from the last step and search in the
    1468 !--          correct direction
    1469              li = li_bnd
    1470              IF ( rib_tab(li) - surf%rib(m) > 0.0_wp )  THEN
    1471                 DO  WHILE ( rib_tab(li-1) - surf%rib(m) > 0.0_wp  .AND.  li > 0 )
    1472                    li = li-1
    1473                 ENDDO
    1474              ELSE
    1475                 DO  WHILE ( rib_tab(li) - surf%rib(m) < 0.0_wp                 &
    1476                            .AND.  li < num_steps-1 )
    1477                    li = li+1
    1478                 ENDDO
    1479              ENDIF
    1480              li_bnd = li
    1481 
    1482 !
    1483 !--          Linear interpolation to find the correct value of z/L
    1484              surf%ol(m) = ( ol_tab(li-1) + ( ol_tab(li) - ol_tab(li-1) )       &
    1485                          / (  rib_tab(li) - rib_tab(li-1) )                    &
    1486                          * ( surf%rib(m) - rib_tab(li-1) ) )
    1487 
    1488           ENDDO
    1489 
    1490        ELSEIF ( TRIM( most_method ) == 'circular' )  THEN
    1491 
    1492           IF ( .NOT. humidity )  THEN
    1493              !$OMP PARALLEL DO PRIVATE( z_mo )
    1494              DO  m = 1, surf%ns
    1495 
    1496                 z_mo = surf%z_mo(m)
    1497 
    1498                 surf%ol(m) =  ( surf%pt1(m) *  surf%us(m)**2 ) /               &
    1499                                        ( kappa * g *                           &
    1500                                          surf%ts(m) + 1E-30_wp )
    1501 !
    1502 !--             Limit the value range of the Obukhov length.
    1503 !--             This is necessary for very small velocities (u,v --> 1), because
    1504 !--             the absolute value of ol can then become very small, which in
    1505 !--             consequence would result in very large shear stresses and very
    1506 !--             small momentum fluxes (both are generally unrealistic).
    1507                 IF ( ( z_mo / ( surf%ol(m) + 1E-30_wp ) ) < zeta_min )         &
    1508                    surf%ol(m) = z_mo / zeta_min
    1509                 IF ( ( z_mo / ( surf%ol(m) + 1E-30_wp ) ) > zeta_max )         &
    1510                    surf%ol(m) = z_mo / zeta_max
    1511 
    1512              ENDDO
    1513           ELSEIF ( bulk_cloud_model  .OR.  cloud_droplets )  THEN
    1514              !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
    1515              DO  m = 1, surf%ns
    1516 
    1517                 i   = surf%i(m)           
    1518                 j   = surf%j(m)
    1519                 k   = surf%k(m)
    1520 
    1521                 z_mo = surf%z_mo(m)
    1522 
    1523 
    1524                 surf%ol(m) =  ( surf%vpt1(m) * surf%us(m)**2 ) /               &
    1525                     ( kappa * g * ( surf%ts(m) +                               &
    1526                          0.61_wp * surf%pt1(m) * surf%us(m)                    &
    1527                        + 0.61_wp * surf%qv1(m) * surf%ts(m) -                  &
    1528                                    surf%ts(m)  * ql(k,j,i) ) + 1E-30_wp )
    1529 !
    1530 !--             Limit the value range of the Obukhov length.
    1531 !--             This is necessary for very small velocities (u,v --> 1), because
    1532 !--             the absolute value of ol can then become very small, which in
    1533 !--             consequence would result in very large shear stresses and very
    1534 !--             small momentum fluxes (both are generally unrealistic).
    1535                 IF ( ( z_mo / ( surf%ol(m) + 1E-30_wp ) ) < zeta_min )         &
    1536                    surf%ol(m) = z_mo / zeta_min
    1537                 IF ( ( z_mo / ( surf%ol(m) + 1E-30_wp ) ) > zeta_max )         &
    1538                    surf%ol(m) = z_mo / zeta_max
    1539 
    1540              ENDDO
    1541           ELSE
    1542 
    1543              !$OMP PARALLEL DO PRIVATE( z_mo )
    1544              DO  m = 1, surf%ns
    1545 
    1546                 z_mo = surf%z_mo(m)
    1547 
    1548                 surf%ol(m) =  ( surf%vpt1(m) *  surf%us(m)**2 ) /              &
    1549                     ( kappa * g * ( surf%ts(m) + 0.61_wp * surf%pt1(m) *       &
    1550                                     surf%qs(m) + 0.61_wp * surf%qv1(m) *       &
    1551                                     surf%ts(m) ) + 1E-30_wp )
    1552 
    1553 !
    1554 !--             Limit the value range of the Obukhov length.
    1555 !--             This is necessary for very small velocities (u,v --> 1), because
    1556 !--             the absolute value of ol can then become very small, which in
    1557 !--             consequence would result in very large shear stresses and very
    1558 !--             small momentum fluxes (both are generally unrealistic).
    1559                 IF ( ( z_mo / ( surf%ol(m) + 1E-30_wp ) ) < zeta_min )         &
    1560                    surf%ol(m) = z_mo / zeta_min
    1561                 IF ( ( z_mo / ( surf%ol(m) + 1E-30_wp ) ) > zeta_max )         &
    1562                    surf%ol(m) = z_mo / zeta_max
    1563 
    1564              ENDDO
    1565 
    1566           ENDIF
    1567 
    1568        ENDIF
     1190       ENDDO
    15691191
    15701192    END SUBROUTINE calc_ol
  • palm/trunk/SOURCE/write_restart_data_mod.f90

    r3655 r3668  
    163163
    164164
    165        binary_version_global = '4.7'
     165       binary_version_global = '4.8'
    166166
    167167       CALL wrd_write_string( 'binary_version_global' )
     
    466466       WRITE ( 14 )  momentum_advec
    467467
    468        CALL wrd_write_string( 'most_method' )
    469        WRITE ( 14 )  most_method
    470 
    471468       CALL wrd_write_string( 'netcdf_precision' )
    472469       WRITE ( 14 )  netcdf_precision
Note: See TracChangeset for help on using the changeset viewer.