Changeset 3668
- Timestamp:
- Jan 14, 2019 12:49:24 PM (6 years ago)
- Location:
- palm/trunk
- Files:
-
- 12 edited
-
SCRIPTS/.csd.config.default (modified) (10 diffs)
-
SCRIPTS/palm_csd (modified) (31 diffs)
-
SCRIPTS/palm_csd_files/palm_csd_canopy_generator.py (modified) (11 diffs)
-
SCRIPTS/palm_csd_files/palm_csd_netcdf_interface.py (modified) (3 diffs)
-
SCRIPTS/palm_csd_files/palm_csd_tools.py (modified) (3 diffs)
-
SOURCE/check_parameters.f90 (modified) (2 diffs)
-
SOURCE/land_surface_model_mod.f90 (modified) (4 diffs)
-
SOURCE/modules.f90 (modified) (2 diffs)
-
SOURCE/parin.f90 (modified) (3 diffs)
-
SOURCE/read_restart_data_mod.f90 (modified) (3 diffs)
-
SOURCE/surface_layer_fluxes_mod.f90 (modified) (9 diffs)
-
SOURCE/write_restart_data_mod.f90 (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SCRIPTS/.csd.config.default
r3567 r3668 25 25 # ----------------- 26 26 # $Id$ 27 # Updated configuration 28 # 29 # 3567 2018-11-27 13:59:21Z maronga 27 30 # Initial revisions 28 #29 #30 #31 #32 31 # 33 32 # Description: … … 52 51 institution = Institute of Meteorology and Climatology, Leibniz University Hannover 53 52 palm_version = 6.0 54 version = 155 53 rotation_angle = 0.0 56 54 references = text … … 58 56 59 57 [settings] 60 filename_out = /ldata2/MOSAIK/showcase_python 61 lai_ season = spring62 lai_high_vegetation_default = 5.063 lai_low_vegetation_default = 1.058 lai_roof_intensive = 1.5 59 lai_roof_extensive = 3.0 60 lai_high_vegetation_default = 6.0 61 lai_low_vegetation_default = 3.0 64 62 lai_alpha = 5.0 65 63 lai_beta = 3.0 66 tree_height_default = 10.064 patch_height_default = 10.0 67 65 bridge_width = 3.0 68 66 debug_mode = False 69 67 season = "summer" 70 68 71 69 [input_01] … … 79 77 file_lon = Berlin_CoordinatesLatLon_x_15m_DLR.nc 80 78 file_zt = Berlin_terrain_height_15m_DLR.nc 81 file_buildings_2d = Berlin_building_height_15m_ LOD2_DLR.nc79 file_buildings_2d = Berlin_building_height_15m_DLR.nc 82 80 file_building_id = Berlin_building_id_15m_DLR.nc 83 81 file_building_type = Berlin_building_type_15m_DLR.nc 84 82 file_bridges_2d = Berlin_bridges_height_15m_DLR.nc 85 83 file_bridges_id = Berlin_bridges_id_15m_DLR.nc 84 file_lai = Berlin_leaf_area_index_15m_DLR_WANG_summer.nc 86 85 file_vegetation_type = Berlin_vegetation_type_15m_DLR.nc 87 86 file_vegetation_height = Berlin_vegetation_patch_height_15m_DLR.nc … … 91 90 file_street_type = Berlin_street_type_15m_DLR.nc 92 91 file_street_crossings = Berlin_street_crossings_15m_DLR.nc 92 file_tree_height = Berlin_trees_height_clean_15m_DLR.nc 93 file_tree_crown_diameter = Berlin_trees_crown_clean_15m_DLR.nc 94 file_tree_trunk_diameter = Berlin_trees_trunk_clean_15m.nc 95 file_tree_type = Berlin_trees_type_15m_DLR.nc 96 file_patch_height = Berlin_vegetation_patch_height_15m_DLR.nc 97 file_vegetation_on_roofs = Berlin_vegetation_on_roofs_15m_DLR.nc 93 98 94 99 [input_02] … … 102 107 file_lon = Berlin_CoordinatesLatLon_x_10m_DLR.nc 103 108 file_zt = Berlin_terrain_height_10m_DLR.nc 104 file_buildings_2d = Berlin_building_height_10m_ LOD2_DLR.nc109 file_buildings_2d = Berlin_building_height_10m_DLR.nc 105 110 file_building_id = Berlin_building_id_10m_DLR.nc 106 111 file_building_type = Berlin_building_type_10m_DLR.nc 107 112 file_bridges_2d = Berlin_bridges_height_10m_DLR.nc 108 113 file_bridges_id = Berlin_bridges_id_10m_DLR.nc 114 file_lai = Berlin_leaf_area_index_10m_DLR_WANG_summer.nc 109 115 file_vegetation_type = Berlin_vegetation_type_10m_DLR.nc 110 116 file_vegetation_height = Berlin_vegetation_patch_height_10m_DLR.nc … … 114 120 file_street_type = Berlin_street_type_10m_DLR.nc 115 121 file_street_crossings = Berlin_street_crossings_10m_DLR.nc 122 file_tree_height = Berlin_trees_height_clean_10m_DLR.nc 123 file_tree_crown_diameter = Berlin_trees_crown_10m_DLR.nc 124 file_tree_trunk_diameter = Berlin_trees_trunk_clean_10m.nc 125 file_tree_type = Berlin_trees_type_10m_DLR.nc 126 file_patch_height = Berlin_vegetation_patch_height_10m_DLR.nc 127 file_vegetation_on_roofs = Berlin_vegetation_on_roofs_10m_DLR.nc 116 128 117 129 [input_03] 130 path = /ldata2/MOSAIK/Berlin_static_driver_data 131 pixel_size = 1.0 118 132 file_x = Berlin_CoordinatesUTM_x_1m_DLR.nc 119 133 file_y = Berlin_CoordinatesUTM_y_1m_DLR.nc … … 122 136 file_lat = Berlin_CoordinatesLatLon_y_1m_DLR.nc 123 137 file_lon = Berlin_CoordinatesLatLon_x_1m_DLR.nc 124 path = /ldata2/MOSAIK/Berlin_static_driver_data125 pixel_size = 1.0126 138 file_zt = Berlin_terrain_height_1m_DLR.nc 127 file_buildings_2d = Berlin_building_height_1m_ LOD2_DLR.nc139 file_buildings_2d = Berlin_building_height_1m_DLR.nc 128 140 file_building_id = Berlin_building_id_1m_DLR.nc 129 141 file_bridges_2d = Berlin_bridges_height_1m_DLR.nc 130 142 file_bridges_id = Berlin_bridges_id_1m_DLR.nc 131 143 file_building_type = Berlin_building_type_1m_DLR.nc 144 file_lai = Berlin_leaf_area_index_1m_DLR_WANG_summer.nc 132 145 file_vegetation_type = Berlin_vegetation_type_1m_DLR.nc 133 146 file_vegetation_height = Berlin_vegetation_patch_height_1m_DLR.nc … … 137 150 file_street_type = Berlin_street_type_1m_DLR.nc 138 151 file_street_crossings = Berlin_street_crossings_1m_DLR.nc 152 file_tree_height = Berlin_trees_height_clean_1m_DLR.nc 153 file_tree_crown_diameter = Berlin_trees_crown_clean_1m_DLR.nc 154 file_tree_trunk_diameter = Berlin_trees_trunk_clean_1m.nc 155 file_tree_type = Berlin_trees_type_1m_DLR.nc 156 file_patch_height = Berlin_vegetation_patch_height_1m_DLR.nc 157 file_vegetation_on_roofs = Berlin_vegetation_on_roofs_1m_DLR.nc 158 159 [input_04] 160 path = /ldata2/MOSAIK/Berlin_static_driver_data 161 pixel_size = 2.0 162 file_x = Berlin_CoordinatesUTM_x_2m_DLR.nc 163 file_y = Berlin_CoordinatesUTM_y_2m_DLR.nc 164 file_x_UTM = Berlin_CoordinatesUTM_x_2m_DLR.nc 165 file_y_UTM = Berlin_CoordinatesUTM_y_2m_DLR.nc 166 file_lat = Berlin_CoordinatesLatLon_y_2m_DLR.nc 167 file_lon = Berlin_CoordinatesLatLon_x_2m_DLR.nc 168 file_zt = Berlin_terrain_height_2m_DLR.nc 169 file_buildings_2d = Berlin_building_height_2m_DLR.nc 170 file_building_id = Berlin_building_id_2m_DLR.nc 171 file_bridges_2d = Berlin_bridges_height_2m_DLR.nc 172 file_bridges_id = Berlin_bridges_id_2m_DLR.nc 173 file_building_type = Berlin_building_type_2m_DLR.nc 174 file_lai = Berlin_leaf_area_index_2m_DLR_WANG_summer.nc 175 file_vegetation_type = Berlin_vegetation_type_2m_DLR.nc 176 file_vegetation_height = Berlin_vegetation_patch_height_2m_DLR.nc 177 file_pavement_type = Berlin_pavement_type_2m_DLR.nc 178 file_water_type = Berlin_water_type_2m_DLR.nc 179 file_soil_type = Berlin_soil_type_2m_DLR.nc 180 file_street_type = Berlin_street_type_2m_DLR.nc 181 file_street_crossings = Berlin_street_crossings_2m_DLR.nc 182 file_tree_height = Berlin_trees_height_clean_2m_DLR.nc 183 file_tree_crown_diameter = Berlin_trees_crown_clean_2m_DLR.nc 184 file_tree_trunk_diameter = Berlin_trees_trunk_clean_2m.nc 185 file_tree_type = Berlin_trees_type_2m_DLR.nc 186 file_patch_height = Berlin_vegetation_patch_height_2m_DLR.nc 187 file_vegetation_on_roofs = Berlin_vegetation_on_roofs_2m_DLR.nc 188 139 189 140 190 [output] 141 191 path = /ldata2/MOSAIK/ 142 file_out = winter_iop1_test192 file_out = showcase_berlin 143 193 version = 1 144 194 145 195 [domain_root] 146 196 pixel_size = 15.0 147 origin_x = 0 .0148 origin_y = 0 .0149 nx = 311 150 ny = 257 197 origin_x = 0 198 origin_y = 0 199 nx = 3119 200 ny = 2573 151 201 buildings_3d = False 152 202 dz = 15.0 153 allow_high_vegetation = True203 allow_high_vegetation = False 154 204 generate_vegetation_patches = True 155 use_palm_z_axis= True205 use_palm_z_axis= False 156 206 interpolate_terrain = False 157 207 domain_parent 208 vegetation_on_roofs = False 209 street_trees = True 158 210 159 211 [domain_N02] … … 170 222 interpolate_terrain = True 171 223 domain_parent = root 224 vegetation_on_roofs = False 225 street_trees = True -
palm/trunk/SCRIPTS/palm_csd
r3629 r3668 25 25 # ----------------- 26 26 # $Id$ 27 # Various improvements and bugfixes 28 # 29 # 3629 2018-12-13 12:18:54Z maronga 27 30 # Added canopy generator calls. Some improvements 28 31 # … … 72 75 73 76 # Definition of settings 74 global settings_filename_out 77 75 78 global settings_bridge_width 76 79 global settings_lai_roof_intensive … … 100 103 global global_site 101 104 global global_source 102 global global_version 105 106 global path_out 107 global filename_out 108 global version_out 103 109 104 110 … … 158 164 global zt_min 159 165 160 settings_filename_out = "default"161 166 settings_bridge_width = 3.0 162 167 settings_lai_roof_intensive = 0.5 … … 185 190 global_site = "" 186 191 global_source = "" 187 global_version = 1 188 192 193 path_out = "" 194 version_out = 1 195 filename_out = "default" 196 189 197 domain_names = [] 190 198 domain_px = [] … … 260 268 global_site = config.get(read_tmp, 'site') 261 269 global_source = config.get(read_tmp, 'source') 262 global_version = float(config.get(read_tmp, 'version')) 270 263 271 264 272 if ( read_tmp == 'settings' ): 265 settings_filename_out = config.get(read_tmp, 'filename_out')266 273 settings_lai_roof_intensive = config.get(read_tmp, 'lai_roof_intensive') 267 274 settings_lai_roof_extensive = config.get(read_tmp, 'lai_roof_extensive') … … 273 280 settings_lai_alpha = float(config.get(read_tmp, 'lai_alpha')) 274 281 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')) 275 287 276 288 if ( read_tmp.split("_")[0] == 'domain' ): … … 322 334 input_file_street_crossings.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_street_crossings')) 323 335 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) 325 342 input_file_tree_height.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_tree_height')) 326 343 input_file_tree_trunk_diameter.append(config.get(read_tmp, 'path') + "/" + config.get(read_tmp, 'file_tree_trunk_diameter')) … … 433 450 # Calculate indices and input files 434 451 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]) 436 453 if domain_parent[i] is not None: 437 454 ii_parent.append(domain_names.index(domain_parent[i])) … … 447 464 # Create NetCDF output file and set global attributes 448 465 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) 450 467 451 468 … … 468 485 del zt_all[:] 469 486 487 print( "Shift terrain height by -" + str(zt_min)) 470 488 for i in range(0,ndomains): 471 489 … … 475 493 y = nc_read_from_file_1d(input_file_y[ii[i]], "y", domain_y0[i], domain_y1[i]) 476 494 477 print( "Shift terrain height by -" + str(zt_min))478 495 zt = zt - zt_min 496 497 nc_write_global_attribute(filename[i], 'origin_z', float(zt_min)) 479 498 480 499 # If necessary, interpolate parent domain terrain height on child domain grid and blend the two … … 494 513 495 514 # Interpolate array and bring to PALM grid of child domain 515 496 516 zt_ip = interpolate_2d(zt_parent,tmp_x,tmp_y,x,y) 497 517 zt_ip = bring_to_palm_grid(zt_ip,x,y,domain_dz[ii_parent[i]]) … … 580 600 581 601 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"] 583 603 building_type = np.where(building_type < 1,defaultvalues["building_type"],building_type) 584 604 … … 594 614 building_type = np.where((building_type == fillvalues["building_type"]) & (buildings_2d != fillvalues["buildings_2d"]),defaultvalues["building_type"],building_type) 595 615 print("Data check #2 " + str(check_arrays_2(buildings_2d,building_type,fillvalues["buildings_2d"],fillvalues["building_type"]))) 616 596 617 597 618 nc_write_to_file_2d(filename[i], 'buildings_2d', buildings_2d, datatypes["buildings_2d"],'y','x',fillvalues["buildings_2d"]) … … 606 627 nc_write_attribute(filename[i], 'building_id', 'long_name', 'building id') 607 628 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]) 609 630 nc_write_attribute(filename[i], 'building_id', 'coordinates', 'E_UTM N_UTM lon lat') 610 631 nc_write_attribute(filename[i], 'building_id', 'grid_mapping', 'E_UTM N_UTM lon lat') … … 616 637 nc_write_attribute(filename[i], 'building_type', 'coordinates', 'E_UTM N_UTM lon lat') 617 638 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))) 618 641 619 642 del buildings_2d … … 705 728 # #6 Remove water for pixels with buildings 706 729 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) 712 737 713 738 # Check for consistency and fill empty fields with default vegetation type … … 717 742 vegetation_type = np.where(consistency_array == 0,defaultvalues["vegetation_type"],vegetation_type) 718 743 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 720 752 # Create surface_fraction array 721 753 x = nc_read_from_file_2d_all(filename[i], 'x') … … 734 766 nc_write_attribute(filename[i], 'surface_fraction', 'res_orig', domain_px[i]) 735 767 del surface_fraction 736 737 738 # Correct vegetation_type when a vegetation height is available and is indicative of low vegeetation739 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 744 768 745 769 nc_write_to_file_2d(filename[i], 'vegetation_type', vegetation_type, datatypes["vegetation_type"],'y','x',fillvalues["vegetation_type"]) … … 779 803 780 804 # 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 782 806 if domain_3d[i]: 783 807 if np.any(building_type != fillvalues["building_type"]): … … 853 877 # Read tree data and create LAD and BAD arrays using the canopy generator 854 878 for 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 907 for i in range(0,ndomains): 855 908 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 859 917 # Save lai data as default for low and high vegetation 860 918 lai_low = lai 861 919 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 882 921 # Read all tree parameters from file 883 922 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 885 932 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]) 886 933 tree_type = nc_read_from_file_2d(input_file_tree_type[ii[i]], 'Band1', domain_x0[i], domain_x1[i], domain_y0[i], domain_y1[i]) … … 889 936 # Remove missing values from the data. Reasonable values will be set by the tree generator 890 937 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)892 938 tree_trunk_diameter = np.where( (tree_trunk_diameter == 0.0) | (tree_trunk_diameter == -1.0) ,fillvalues["tree_data"],tree_trunk_diameter) 893 939 tree_type = np.where( (tree_type == 0.0) | (tree_type == -1.0) ,fillvalues["tree_data"],tree_type) … … 907 953 # For zero-height patches, set vegetation_type to short grass and remove these pixels from the patch height array 908 954 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 911 957 912 958 max_tree_height = max(tree_height.flatten()) … … 914 960 915 961 if ( (max_tree_height != fillvalues["tree_data"]) | (max_patch_height == fillvalues["tree_data"]) ): 916 962 917 963 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"]) 918 964 … … 925 971 bad[k,:,:] = np.where(buildings_2d == fillvalues["buildings_2d"],bad[k,:,:],fillvalues["tree_data"]) 926 972 tree_ids[k,:,:] = np.where(buildings_2d == fillvalues["buildings_2d"],tree_ids[k,:,:],fillvalues["tree_data"]) 927 973 974 del buildings_2d 928 975 929 976 nc_write_dimension(filename[i], 'zlad', zlad, datatypes["tree_data"]) … … 949 996 nc_write_attribute(filename[i], 'tree_id', 'grid_mapping', 'E_UTM N_UTM lon lat') 950 997 951 del buildings_2d,lai, lad, bad, tree_ids, zlad952 953 del vegetation_pars, tree_height, tree_crown_diameter, tree_trunk_diameter, tree_type, patch_height, vegetation_type,x, y998 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 954 1001 955 1002 … … 961 1008 vegetation_type = nc_read_from_file_2d_all(filename[i], 'vegetation_type') 962 1009 lad = nc_read_from_file_3d_all(filename[i], 'lad') 1010 zlad = nc_read_from_file_1d_all(filename[i], 'zlad') 963 1011 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]) 964 1012 vegetation_pars = nc_read_from_file_3d_all(filename[i], 'vegetation_pars') 965 1013 lai = vegetation_pars[1,:,:] 966 1014 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 value969 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 where ever it is missing, but where a high vegetation LAI was set1015 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 975 1023 patch_height = np.where ( (lai_high != fillvalues["tree_data"]) & (patch_height == fillvalues["tree_data"]), settings_patch_height_default, patch_height) 976 1024 … … 978 1026 patch_height = np.where ( (lad[0,:,:] != fillvalues["tree_data"]), fillvalues["tree_data"], patch_height) 979 1027 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 980 1031 # 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 985 1035 if ( max(patch_height.flatten()) >= (2.0 * domain_dz[i]) ): 986 1036 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) 988 1038 989 1039 lad[0:patch_nz+1,:,:] = np.where( (lad[0:patch_nz+1,:,:] == fillvalues["tree_data"]),lad_patch[0:patch_nz+1,:,:], lad[0:patch_nz+1,:,:] ) … … 991 1041 # Remove high vegetation wherever it is replaced by a leaf area density. This should effectively remove all high vegetation pixels 992 1042 vegetation_type = np.where((lad[0,:,:] != fillvalues["tree_data"]) & (vegetation_type != fillvalues["vegetation_type"]),3,vegetation_type) 993 994 1043 995 1044 # If desired, remove all high vegetation. TODO: check if this is still necessary 996 1045 if domain_high_vegetation[i]: … … 998 1047 999 1048 1000 # Remove LAI from pixels where a LAD was set1001 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) 1002 1051 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 1007 1058 # Overwrite lai in vegetation_parameters 1008 1009 vegetation_pars[1,:,:] = lai_low 1059 vegetation_pars[1,:,:] = np.copy(lai_low) 1010 1060 nc_overwrite_to_file_3d(filename[i], 'vegetation_pars', vegetation_pars) 1011 1061 1012 1062 # Overwrite lad array 1013 1063 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 1071 for 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 25 25 # ----------------- 26 26 # $Id$ 27 # Various improvements and bugfixes 28 # 29 # 3629 2018-12-13 12:18:54Z maronga 27 30 # Initial revision 28 31 # … … 72 75 if valid_pixels[j,i]: 73 76 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="") 76 78 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) 77 79 if ( np.any(lad_loc) != fill ): … … 101 103 102 104 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) 108 110 109 111 del lad_loc, x_loc, y_loc, z_loc, status … … 131 133 default_trees = [] 132 134 #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)) 200 202 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)) 220 222 221 223 # Define fill values … … 237 239 if (season == "summer"): 238 240 tree_lai = default_trees[tree_type].lai_summer 241 239 242 else: 240 243 tree_lai = default_trees[tree_type].lai_winter … … 279 282 280 283 lad_max = tree_lai / (lad_max_part_1[0] + lad_max_part_2[0]) 281 282 284 283 285 … … 298 300 z = np.arange(0,nz*dz,dz) 299 301 z[1:] = z[1:] - 0.5 * dz 300 301 302 302 303 # Define center of the tree position inside the local LAD domain … … 468 469 469 470 470 def process_patch(dz,patch_height, patch_lai,alpha,beta):471 def process_patch(dz,patch_height,max_height_lad,patch_lai,alpha,beta): 471 472 472 473 # Define fill values … … 481 482 pch_index = np.where( (pch_index[:,:] == -1) ,0,pch_index[:,:]) 482 483 483 max_canopy_height = max( patch_height.flatten())484 max_canopy_height = max(max(patch_height.flatten()),max_height_lad) 484 485 485 486 z = np.arange(0,math.floor(max_canopy_height/dz)*dz+2*dz,dz) … … 500 501 int_bpdf = 0.0 501 502 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]): 503 504 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] ) ) 504 505 505 for k in range( 0,pch_index[j,i]):506 for k in range(1,pch_index[j,i]): 506 507 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] 507 508 … … 538 539 # dbh: default trunk diameter at breast height (1.4 m) (m) 539 540 # 540 #541 541 class tree: 542 542 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 25 25 # ----------------- 26 26 # $Id$ 27 # Some improvements and new routines 28 # 29 # 3629 2018-12-13 12:18:54Z maronga 27 30 # Some new routines 28 31 # … … 39 42 40 43 from netCDF4 import Dataset 44 45 def 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 41 63 42 64 def nc_read_from_file_2d(filename, varname, x0, x1, y0, y1): … … 174 196 return 0 175 197 198 def 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 176 211 177 212 def nc_write_dimension(filename,varname,array,datatype): -
palm/trunk/SCRIPTS/palm_csd_files/palm_csd_tools.py
r3567 r3668 25 25 # ----------------- 26 26 # $Id$ 27 # Minor changes 28 # 29 # 3567 2018-11-27 13:59:21Z maronga 27 30 # Initial revisions 28 #29 #30 #31 31 # 32 32 # … … 77 77 78 78 def 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)) 82 81 83 82 return array_ip 84 83 85 84 … … 136 135 return result 137 136 137 def 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 138 151 139 152 def check_consistency_4(array1,array2,array3,array4,fill1,fill2,fill3,fill4): -
palm/trunk/SOURCE/check_parameters.f90
r3655 r3668 25 25 ! ----------------- 26 26 ! $Id$ 27 ! most_method removed 28 ! 29 ! 3655 2019-01-07 16:51:22Z knoop 27 30 ! Formatting 28 31 ! … … 3790 3793 ENDIF 3791 3794 3792 !3793 !-- Check for valid setting of most_method3794 IF ( TRIM( most_method ) /= 'circular' .AND. &3795 TRIM( most_method ) /= 'newton' .AND. &3796 TRIM( most_method ) /= 'lookup' ) THEN3797 message_string = 'most_method = "' // TRIM( most_method ) // &3798 '" is unknown'3799 CALL message( 'check_parameters', 'PA0416', 1, 2, 0, 6, 0 )3800 ENDIF3801 3795 3802 3796 ! -
palm/trunk/SOURCE/land_surface_model_mod.f90
r3655 r3668 25 25 ! ----------------- 26 26 ! $Id$ 27 ! Removed most_method 28 ! 29 ! 3655 2019-01-07 16:51:22Z knoop 27 30 ! nopointer option removed 28 31 ! … … 1324 1327 1325 1328 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 1328 1330 1329 1331 … … 1446 1448 1447 1449 IF ( TRIM( surface_type ) == 'water' ) THEN 1448 1449 IF ( TRIM( most_method ) == 'lookup' ) THEN1450 WRITE( message_string, * ) 'surface_type = ', surface_type, &1451 ' is not allowed in combination with ', &1452 'most_method = ', most_method1453 CALL message( 'lsm_check_parameters', 'PA0414', 1, 2, 0, 6, 0 )1454 ENDIF1455 1450 1456 1451 IF ( water_type == 0 ) THEN … … 1531 1526 1532 1527 IF ( TRIM( surface_type ) == 'netcdf' ) THEN 1533 IF ( ANY( water_type_f%var /= water_type_f%fill ) .AND. &1534 TRIM( most_method ) == 'lookup' ) THEN1535 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 ENDIF1540 1528 ! 1541 1529 !-- MS: Some problme here, after calling message everythings stucks at -
palm/trunk/SOURCE/modules.f90
r3648 r3668 25 25 ! ----------------- 26 26 ! $Id$ 27 ! Removed most_method 28 ! 29 ! 3648 2019-01-02 16:35:46Z suehring 27 30 ! -surface_data_output +surface_output 28 31 ! … … 1090 1093 CHARACTER (LEN=8) :: coupling_char = '' !< appended to filenames in coupled or nested runs ('_O': ocean PE, 1091 1094 !< '_NV': vertically nested atmosphere PE, '_N##': PE of nested domain ## 1092 CHARACTER (LEN=8) :: most_method = 'newton' !< namelist parameter1093 1095 CHARACTER (LEN=10) :: run_date = ' ' !< date of simulation run 1094 1096 CHARACTER (LEN=8) :: run_time = ' ' !< time of simulation run -
palm/trunk/SOURCE/parin.f90
r3649 r3668 25 25 ! ----------------- 26 26 ! $Id$ 27 ! Removed most_method 28 ! 29 ! 3649 2019-01-02 16:52:21Z suehring 27 30 ! Delete debug-print statements 28 31 ! … … 573 576 loop_optimization, lsf_exception, masking_method, mg_cycles, & 574 577 mg_switch_to_pe0_level, mixing_length_1d, momentum_advec, & 575 most_method, &576 578 netcdf_precision, neutral, ngsrb, & 577 579 nsor, nsor_ini, nudging, nx, ny, nz, ocean_mode, omega, & … … 644 646 loop_optimization, lsf_exception, masking_method, mg_cycles, & 645 647 mg_switch_to_pe0_level, mixing_length_1d, momentum_advec, & 646 most_method, &647 648 netcdf_precision, neutral, ngsrb, & 648 649 nsor, nsor_ini, nudging, nx, ny, nz, ocean_mode, omega, & -
palm/trunk/SOURCE/read_restart_data_mod.f90
r3655 r3668 25 25 ! ----------------- 26 26 ! $Id$ 27 ! Removed most_method and increased binary version 28 ! 29 ! 3655 2019-01-07 16:51:22Z knoop 27 30 ! Implementation of the PALM module interface 28 31 ! … … 203 206 READ ( 13 ) version_on_file 204 207 205 binary_version_global = '4. 7'208 binary_version_global = '4.8' 206 209 IF ( TRIM( version_on_file ) /= TRIM( binary_version_global ) ) THEN 207 210 WRITE( message_string, * ) 'version mismatch concerning ', & … … 488 491 CASE ( 'momentum_advec' ) 489 492 READ ( 13 ) momentum_advec 490 CASE ( 'most_method' )491 READ ( 13 ) most_method492 493 CASE ( 'netcdf_precision' ) 493 494 READ ( 13 ) netcdf_precision -
palm/trunk/SOURCE/surface_layer_fluxes_mod.f90
r3655 r3668 26 26 ! ----------------- 27 27 ! $Id$ 28 ! Removed methods "circular" and "lookup" 29 ! 30 ! 3655 2019-01-07 16:51:22Z knoop 28 31 ! OpenACC port for SPEC 29 32 ! … … 228 231 ! ------------ 229 232 !> 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 237 234 !> 238 235 !> @todo (re)move large_scale_forcing actions … … 261 258 intermediate_timestep_count, intermediate_timestep_count_max, & 262 259 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, & 264 261 run_coupled, surface_pressure, simulated_time, terminate_run, & 265 262 time_since_reference_point, urban_surface, & … … 400 397 !-- First, calculate the new Obukhov length, then new friction velocity, 401 398 !-- 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. 409 400 !-- First call for default-type horizontal surfaces, for natural- and 410 401 !-- 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 503 416 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 504 447 ! 505 448 !-- Treat downward-facing horizontal surfaces. Note, so far, these are … … 549 492 !-- flux in land-surface model in case of aero_resist_kray is not true. 550 493 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 594 514 ! 595 515 !-- No ts is required, so scaling parameters and Obukhov length do not need … … 652 572 !-- flux in land-surface model in case of aero_resist_kray is not true. 653 573 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 697 594 ELSE 698 595 DO l = 2, 3 … … 841 738 IF ( surf_lsm_h%ns >= 1 ) surf_lsm_h%ol = 1.0E10_wp 842 739 IF ( surf_usm_h%ns >= 1 ) surf_usm_h%ol = 1.0E10_wp 843 ENDIF844 845 IF ( most_method == 'lookup' ) THEN846 847 !848 !-- Check for roughness heterogeneity. In that case terminate run and849 !-- inform user. Check for both, natural and non-natural walls.850 IF ( surf_def_h(0)%ns >= 1 ) THEN851 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 ) ) THEN853 terminate_run_l = .TRUE.854 ENDIF855 ENDIF856 IF ( surf_lsm_h%ns >= 1 ) THEN857 IF ( MINVAL( surf_lsm_h%z0h ) /= MAXVAL( surf_lsm_h%z0h ) .OR. &858 MINVAL( surf_lsm_h%z0 ) /= MAXVAL( surf_lsm_h%z0 ) ) THEN859 terminate_run_l = .TRUE.860 ENDIF861 ENDIF862 IF ( surf_usm_h%ns >= 1 ) THEN863 IF ( MINVAL( surf_usm_h%z0h ) /= MAXVAL( surf_usm_h%z0h ) .OR. &864 MINVAL( surf_usm_h%z0 ) /= MAXVAL( surf_usm_h%z0 ) ) THEN865 terminate_run_l = .TRUE.866 ENDIF867 ENDIF868 !869 !-- Check roughness homogeneity between differt surface types.870 IF ( surf_lsm_h%ns >= 1 .AND. surf_def_h(0)%ns >= 1 ) THEN871 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 ) ) THEN873 terminate_run_l = .TRUE.874 ENDIF875 ENDIF876 IF ( surf_usm_h%ns >= 1 .AND. surf_def_h(0)%ns >= 1 ) THEN877 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 ) ) THEN879 terminate_run_l = .TRUE.880 ENDIF881 ENDIF882 IF ( surf_usm_h%ns >= 1 .AND. surf_lsm_h%ns >= 1 ) THEN883 IF ( MINVAL( surf_usm_h%z0h ) /= MAXVAL( surf_lsm_h%z0h ) .OR. &884 MINVAL( surf_usm_h%z0 ) /= MAXVAL( surf_lsm_h%z0 ) ) THEN885 terminate_run_l = .TRUE.886 ENDIF887 ENDIF888 889 890 #if defined( __parallel )891 !892 !-- Make a logical OR for all processes. Force termiation of model if result893 !-- is TRUE894 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 #else898 terminate_run = terminate_run_l899 #endif900 901 IF ( terminate_run ) THEN902 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 ENDIF907 908 ALLOCATE( zeta_tmp(0:num_steps-1) )909 910 !911 !-- Use the lowest possible value for z_mo912 k = nzb913 z_mo = zu(k+1) - zw(k)914 915 !916 !-- Calculate z/L range from zeta_stretch to zeta_max using 90% of the917 !-- available steps (num_steps). The calculation is done with negative918 !-- values of zeta in order to simplify the stretching in the free919 !-- convection limit for the remaining 10% of steps.920 zeta_tmp(0) = - zeta_max921 num_steps_n = ( num_steps * 9 / 10 ) - 1922 zeta_step = (zeta_max - zeta_stretch) / REAL(num_steps_n)923 924 DO li = 1, num_steps_n925 zeta_tmp(li) = zeta_tmp(li-1) + zeta_step926 ENDDO927 928 !929 !-- Calculate stretching factor for the free convection range930 DO WHILE ( ABS( (regr-regr_old) / regr_old ) > 1.0E-10_wp )931 regr_old = regr932 regr = ( 1.0_wp - ( -zeta_min / zeta_step ) * ( 1.0_wp - regr ) &933 )**( 10.0_wp / REAL(num_steps) )934 ENDDO935 936 !937 !-- Calculate z/L range from zeta_min to zeta_stretch938 DO li = num_steps_n+1, num_steps-1939 zeta_tmp(li) = zeta_tmp(li-1) + zeta_step940 zeta_step = zeta_step * regr941 ENDDO942 943 !944 !-- Save roughness lengths to temporary variables945 IF ( surf_def_h(0)%ns >= 1 ) THEN946 z0h_min = surf_def_h(0)%z0h(1)947 z0_min = surf_def_h(0)%z0(1)948 ELSEIF ( surf_lsm_h%ns >= 1 ) THEN949 z0h_min = surf_lsm_h%z0h(1)950 z0_min = surf_lsm_h%z0(1)951 ELSEIF ( surf_usm_h%ns >= 1 ) THEN952 z0h_min = surf_usm_h%z0h(1)953 z0_min = surf_usm_h%z0(1)954 ENDIF955 !956 !-- Calculate lookup table for the Richardson number versus Obukhov length957 !-- The Richardson number (rib) is defined depending on the choice of958 !-- boundary conditions for temperature959 IF ( ibc_pt_b == 1 ) THEN960 DO li = 0, num_steps-1961 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 )**3966 ENDDO967 ELSE968 DO li = 0, num_steps-1969 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 )**2978 ENDDO979 ENDIF980 981 !982 !-- Determine minimum values of rib in the lookup table. Set upper limit983 !-- 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 )988 740 ENDIF 989 741 … … 1256 1008 ol_u !< Upper bound of L for Newton iteration 1257 1009 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 1286 1025 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 & 1296 1046 * surf%qv1(m) ) * surf%shf(m) + 0.61_wp & 1297 1047 * surf%pt1(m) * surf%qsws(m) ) * & … … 1299 1049 ( surf%uvw_abs(m)**3 * surf%vpt1(m) * kappa**2 & 1300 1050 + 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 ) 1302 1145 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 ) 1317 1165 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 1341 1169 1342 1170 ! 1343 1171 !-- 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 1349 1180 ENDIF 1350 1181 ! 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 1449 1189 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 1569 1191 1570 1192 END SUBROUTINE calc_ol -
palm/trunk/SOURCE/write_restart_data_mod.f90
r3655 r3668 163 163 164 164 165 binary_version_global = '4. 7'165 binary_version_global = '4.8' 166 166 167 167 CALL wrd_write_string( 'binary_version_global' ) … … 466 466 WRITE ( 14 ) momentum_advec 467 467 468 CALL wrd_write_string( 'most_method' )469 WRITE ( 14 ) most_method470 471 468 CALL wrd_write_string( 'netcdf_precision' ) 472 469 WRITE ( 14 ) netcdf_precision
Note: See TracChangeset
for help on using the changeset viewer.
