Changeset 3727
 Timestamp:
 Feb 8, 2019 2:52:10 PM (4 years ago)
 Location:
 palm/trunk
 Files:

 4 edited
Legend:
 Unmodified
 Added
 Removed

palm/trunk/SCRIPTS/.palm.iofiles
r3707 r3727 59 59 SURFACE_DATA_BIN* out:tr * $base_data/$run_identifier/OUTPUT _surf bin 60 60 SURFACE_DATA_AV_BIN* out:tr * $base_data/$run_identifier/OUTPUT _av_surf bin 61 SURFACE_DATA_NETCDF* out:tr * $base_data/$run_identifier/OUTPUT _surf nc 62 SURFACE_DATA_AV_NETCDF* out:tr * $base_data/$run_identifier/OUTPUT _av_surf nc 61 63 # 62 64 VIRTUAL_MEAS_BIN* out * $base_data/$run_identifier/OUTPUT _vmeas bin 
palm/trunk/SOURCE/Makefile
r3700 r3727 25 25 #  26 26 # $Id$ 27 # surface_data_output_mod depends on netcdf_interface_mod 28 # 29 # 3700 20190126 17:03:42Z knoop 27 30 # Rename surface_output_mod into surface_data_output_mod 28 # 31 # 29 32 # 3637 20181220 01:51:36Z knoop 30 33 # Implementation of the PALM module interface 31 # 34 # 32 35 # 3634 20181218 12:31:28Z knoop 33 36 # OpenACC port for SPEC 34 # 37 # 35 38 # 3579 20181129 15:32:39Z suehring 36 39 # Dependency for check_parameters on nesting_offl_mod added 37 # 40 # 38 41 # 3569 20181127 17:03:40Z kanani 39 42 # dom_dwd_user, Schrempf: … … 43 46 # 3525 20181114 16:06:14Z kanani 44 47 # Changes related to cleanup of biometeorology (dom_dwd_user) 45 # 48 # 46 49 # 3522 20181113 12:14:36Z suehring 47 50 # Dependencies for virtual measurement module added 48 # 51 # 49 52 # 3494 20181106 14:51:27Z suehring 50 53 # Surface output revised 51 # 54 # 52 55 # 3474 20181030 21:07:39Z kanani 53 56 # Add virtual measurement module 54 # 57 # 55 58 # 3472 20181030 20:43:50Z suehring 56 59 # Add indoor model (kanani, srissman, tlang), 57 60 # minor formatting 58 # 61 # 59 62 # 3467 20181030 19:05:21Z suehring 60 63 # Implementation of a new aerosol module salsa. 61 # 62 # 64 # 65 # 63 66 # 3458 20181030 14:51:23Z kanani 64 67 # from chemistry branch r3443, banzhafs, Russo, forkel, basit: … … 67 70 # Added chemistry emission module 68 71 # chemistry_model_mod added to flow_statistics 69 # 72 # 70 73 # 3448 20181029 18:14:31Z kanani 71 74 # Adjustment of biometeorology dependencies 72 # 75 # 73 76 # 3436 20181026 18:35:15Z gronemeier 74 77 # Add surface_mod to user_data_output_mask 75 # 78 # 76 79 # 3435 20181026 18:25:44Z gronemeier 77 80 #  Add surface_mod to data_output_mask 78 81 #  Add chemistry_model_mod and surface_mod to init_masks 79 # 82 # 80 83 # 3421 20181024 18:39:32Z gronemeier 81 84 # Add netcdf_data_input_mod to netcdf_interface_mod 82 85 # bugfix: add dependencies to chemistry_model_mod 83 86 # Add module for surface data output 84 # 87 # 85 88 # 3381 20181019 13:09:06Z raasch 86 89 # dependencies for ocean_mod fixed 87 # 90 # 88 91 # 3355 20181016 14:03:34Z knoop 89 92 # Add module for offline nesting; … … 91 94 # Bugfix, missing dependency for turbulence generator in init_3d_model; 92 95 # Some formatting ajdustments 93 # 96 # 94 97 # 3343 20181015 10:38:52Z suehring 95 98 # (from branch resler) 96 99 # Add biometeorology 97 # 98 # 100 # 101 # 99 102 # 3322 20181009 10:02:39Z kanani 100 103 # Formatting and cleanup 101 # 104 # 102 105 # 3298 20181002 12:21:11Z kanani 103 106 # Added missing dependencies and replaced blanks with tabs (forkel) … … 108 111 # changes related to modularization of the ocean mode, 109 112 # bugfix: dependency to advec_ws was missed in chemistry_model_mod 110 # 113 # 111 114 # 3274 20180924 15:42:55Z knoop 112 115 # Added palm dependency of multi_agent_system_mod, because of mas_last_actions 113 116 # call at the end of palm run 114 # 117 # 115 118 # 3167 20180724 18:17:30Z suehring 116 119 # Bugfix, add missing dependencies for multiagent system 117 # 120 # 118 121 # 3159 20180720 11:20:01Z sward 119 122 # Added multi agent system 120 # 123 # 121 124 # 3130 20180716 11:08:55Z gronemeier 122 125 # add surface_layer_fluxes_mod to turbulence_closure_mod 123 # 126 # 124 127 # 3129 20180716 07:45:13Z gronemeier 125 128 # add turbulence_closure_mod to parin 126 # 129 # 127 130 # 2963 20180412 14:47:44Z suehring 128 # Introduce index for vegetation/wall, pavement/greenwall and water/window 131 # Introduce index for vegetation/wall, pavement/greenwall and water/window 129 132 # surfaces, for clearer access of surface fraction, albedo, emissivity, etc. . 130 # 133 # 131 134 # 2955 20180409 15:14:01Z suehring 132 135 # Add logpoints to measure CPU time of NetCDF data input. 133 # 136 # 134 137 # 2938 20180327 15:52:42Z suehring 135 138 # No initialization of child domains via dynamic input file, except for soil 136 139 # moisture and temperature 137 140 # Apply turbulence generator at noncyclic lateral boundary in nesting case 138 # 141 # 139 142 # 2936 20180327 14:49:27Z suehring 140 143 # Added dependencies for parent and child synchronization 141 # 144 # 142 145 # 2921 20180322 15:05:23Z Giersch 143 146 # date_and_time_mod dependency has been added to read/write_restart_data_mod 144 # 147 # 145 148 # 2918 20180321 15:52:14Z gronemeier 146 # read/write_3d_binary and read/write_var_list has been removed, 149 # read/write_3d_binary and read/write_var_list has been removed, 147 150 # read/write_restart_data_mod, wrd_write_string and 148 # user_read/write_restart_data_mod has been added, dependencies with respect to 151 # user_read/write_restart_data_mod has been added, dependencies with respect to 149 152 # the aforementioned routines have been added/removed 150 153 # 151 154 # 2847 20180302 21:45:58Z suehring 152 155 # Changed format and enforced sorting 153 # 156 # 154 157 # 2817 20180219 16:32:21Z knoop 155 158 # Preliminary gust module interface implemented 156 # 159 # 157 160 # 2802 20180214 16:21:39Z thiele 158 161 # Changed lpm from subroutine to module. 159 162 # Introduce particle transfer in nested models. 160 # 163 # 161 164 # 2773 20180130 14:12:54Z suehring 162 165 # Nesting of chemical species 163 # 166 # 164 167 # 2718 20180102 08:49:38Z maronga 165 168 # Corrected "Former revisions" section 166 # 169 # 167 170 # 2697 20171214 17:57:20Z kanani 168 171 # Bugfix, missing dependencies 169 # 172 # 170 173 # 2696 20171214 17:12:51Z kanani 171 174 # Change in file header (GPL part) 172 175 # Implementation of uv exposure model (FK) 173 # Bugfix, removed loop dependcy for vertical_nesting_mod and 176 # Bugfix, removed loop dependcy for vertical_nesting_mod and 174 177 # turbulence_closure_mod, added depencies for vertical_nesting_mod (TG) 175 178 # implemented turbulence_closure_mod (TG) … … 179 182 # For LSM, add dependency on calc_mean_profile (??) 180 183 # poismg_noopt modularized and renamed into poismg_noopt_mod 181 # add dependencies for netcdf_data_input_mod, calc_mean_profile, 184 # add dependencies for netcdf_data_input_mod, calc_mean_profile, 182 185 # radiation_module_mod, land_surface_model_mod (MS) 183 # 186 # 184 187 # 2608 20171113 14:04:26Z schwenkel 185 # Added diagnostic_quantities_mod 186 # 188 # Added diagnostic_quantities_mod 189 # 187 190 # 2600 20171101 14:11:20Z raasch 188 191 # comment line concerning bound checks removed 189 # 192 # 190 193 # 2599 20171101 13:18:45Z hellstea 191 194 # virtual_flight_mod, synthetic_turbulence_generator_mod and … … 194 197 # 2563 20171019 15:36:10Z Giersch 195 198 # wind_turbine_model_mod and synthetic_turbulence_generator_mod were added to 196 # write_var_list and virtual_flight_mod was deleted from read_var_list 199 # write_var_list and virtual_flight_mod was deleted from read_var_list 197 200 # 198 201 # 2544 20171013 18:09:32Z maronga 199 # Added date_and_time_mod 200 # 202 # Added date_and_time_mod 203 # 201 204 # 2371 20170824 13:01:17Z kanani 202 205 # Corrected dependencies for vertical_nesting_mod 203 # 206 # 204 207 # 2370 20170823 06:11:43Z raasch 205 208 # dependency bugfix for synthetic_turbulence_generator 206 # 209 # 207 210 # 2365 20170821 14:59:59Z kanani 208 211 # Added dependencies for vertical_nesting_mod 209 # 212 # 210 213 # 2339 20170807 13:55:26Z gronemeier 211 214 # corrected timestamp in header 212 # 215 # 213 216 # 2338 20170807 12:15:38Z gronemeier 214 217 # Modularize 1D model 215 # 218 # 216 219 # 2320 20170721 12:47:43Z suehring 217 220 # ls_forcing nudging 218 221 # +large_scale_forcing_nudging 219 # 222 # 220 223 # 2318 20170720 17:27:44Z suehring 221 224 # Add further dependencies on surface_mod 222 # 225 # 223 226 # 2317 20170720 17:27:19Z suehring 224 227 # Added time_integration_spinup 225 # 228 # 226 229 # 2269 20170609 11:57:32Z suehring 227 230 # Add dependency in read_3d_binary 228 # 231 # 229 232 # 2263 20170608 14:59:01Z schwenkel 230 233 # Implemented splitting and merging algorithm 231 # 234 # 232 235 # 2259 20170608 09:09:11Z gronemeier 233 236 # Implemented synthetic turbulence generator … … 235 238 # 2256 20170607 13:58:08Z suehring 236 239 # Remove ring dependency in init_pegrid 237 # 240 # 238 241 # 2238 20170531 16:49:16Z suehring 239 242 # Bugfix, further missing dependency on surface_mod 240 # 243 # 241 244 # 2237 20170531 10:34:53Z suehring 242 # Bugfix, add dependencies on surface_mod for surface_coupler, 243 # plant_canopy_model_mod and ls_forcing_mod 244 # 245 # Bugfix, add dependencies on surface_mod for surface_coupler, 246 # plant_canopy_model_mod and ls_forcing_mod 247 # 245 248 # 2233 20170530 18:08:54Z suehring 246 249 # … … 248 251 # +dependencies for surface_mod 249 252 # wall_fluxes 250 # 253 # 251 254 # 2130 20170124 16:25:39Z raasch 252 255 # dependency for timestep updated … … 254 257 # 2118 20170117 16:38:49Z raasch 255 258 # cuda_fft_interfaces_mod 256 # 259 # 257 260 # 2050 20161108 15:00:55Z gronemeier 258 261 # Implement turbulent outflow method 259 # 262 # 260 263 # 2007 20160824 15:47:17Z kanani 261 264 # urban surface module added, 262 265 # cleaned up some lines (compiler flags/options), which were accidentally 263 266 # added in rev1938 264 # 267 # 265 268 # 1998 20160820 18:45:34Z knoop 266 269 # Bugfix: added netcdf_interface to dependency list for user_init_land_surface … … 268 271 # 1986 20160810 14:07:17Z gronemeier 269 272 # POSIXcalls module added 270 # 273 # 271 274 # 1972 20160726 07:52:02Z maronga 272 275 # Removed some dependencies due to further modularization of land surface model 273 # 276 # 274 277 # 1957 20160707 10:43:48Z suehring 275 278 # flight module added … … 279 282 # 280 283 # 1938 20160613 15:26:05Z hellstea 281 # Some dependency errors corrected 282 # 284 # Some dependency errors corrected 285 # 283 286 # 1934 20160613 09:46:57Z hellstea 284 287 # poismg renamed poismg_noopt, poismg_fast_mod renamed poismg_mod 285 # 288 # 286 289 # 1914 20160526 14:44:07Z witha 287 290 # Added wind_turbine_model_mod … … 294 297 # 295 298 # 1850 20160408 13:29:27Z maronga 296 # Adapted for modularization of microphysics 299 # Adapted for modularization of microphysics 297 300 # Several files renamed > _mod 298 301 # Bugfix for previous commit … … 309 312 # 1833 20160407 14:23:03Z raasch 310 313 # spectrum renamed spectra_mod, depencies for spectra changed 311 # 314 # 312 315 # 1826 20160407 12:01:39Z maronga 313 316 # Renamed radiation_model to radiation_model_mod. 314 317 # Renamed plant_canopy_model to plant_canopy_model_mod. 315 # 318 # 316 319 # 1822 20160407 07:49:42Z hoffmann 317 320 # Tails removed. lpm_release_set removed. calc_precipitation, impact_of_latent_heat 318 # removed. 321 # removed. 319 322 # 320 323 # 1817 20160406 15:44:20Z maronga 321 324 # Renamed land_surface_model to land_surface_model_mod. 322 # 325 # 323 326 # 1808 20160405 19:44:00Z raasch 324 327 # local_flush, local_getenv … … 346 349 # 347 350 # 1762 20160225 12:31:13Z hellstea 348 # +pmc_interface, +pmc routines 351 # +pmc_interface, +pmc routines 349 352 # 350 353 # 1747 20160208 12:25:53Z raasch … … 352 355 # 353 356 # 1691 20151026 16:17:44Z maronga 354 # Replaced prandtl_fluxes with surface_layer_fluxes. Added radiation_model to 357 # Replaced prandtl_fluxes with surface_layer_fluxes. Added radiation_model to 355 358 # prognostic_equations 356 # 359 # 357 360 # 1585 20150430 07:05:52Z maronga 358 361 # Added user_init_radiation.f90 359 # 362 # 360 363 # 1575 20150327 09:56:27Z raasch 361 364 # +poismg_fast … … 363 366 # 1551 20150303 14:18:16Z maronga 364 367 # Bugfix: further adjustments for the land surface model and radiation model 365 # 368 # 366 369 # 1517 20150107 19:12:25Z hoffmann 367 370 # advec_s_bc added to prognostic_equations … … 369 372 # 1500 20141203 17:42:41Z maronga 370 373 # Bugfix: missing adjustments for land surface model and radiation model 371 # 374 # 372 375 # 1496 20141202 17:25:50Z maronga 373 # Added land surface model and radiation model files: land_surface_model, 376 # Added land surface model and radiation model files: land_surface_model, 374 377 # radiation_model, user_init_land_surface 375 # 378 # 376 379 # 1484 20141021 10:53:05Z kanani 377 380 # plant_canopy_modeldependency added for check_parameters, header, init_3d_model, 378 381 # package_parin, read_var_list, user_init_plant_canopy, write_var_list 379 # 382 # 380 383 # 1444 20140802 20:10:32Z letzel 381 384 # bugfix: cpulog added to lpm_advec 382 # 385 # 383 386 # 1404 20140514 09:01:39Z keck 384 387 # bugfix: dependencies added for progress_bar … … 386 389 # 1402 20140509 14:25:13Z raasch 387 390 # progress_bar added 388 # 391 # 389 392 # 1400 20140509 14:03:54Z knoop 390 393 # Added new module random_generator_parallel 391 # 394 # 392 395 # 1380 20140428 12:40:45Z heinze 393 396 # bugfix: mod_particle_attributes added to check_open 394 # nudging added to time_integration 397 # nudging added to time_integration 395 398 # 396 399 # 1374 20140425 12:55:07Z raasch … … 400 403 # Added new module calc_mean_profile, previously in module buoyancy, 401 404 # removed buoyancy dependency from nudging 402 # 405 # 403 406 # 1363 20140417 12:28:49Z keck 404 407 # bugfix: cpulog added to lpm_pack_arrays … … 406 409 # 1361 20140416 15:17:48Z hoffmann 407 410 # cpulog added to microphysics 408 # 411 # 409 412 # 1359 20140411 17:15:14Z hoffmann 410 # mod_particle_attributes added, lpm_sort_arrays removed, 413 # mod_particle_attributes added, lpm_sort_arrays removed, 411 414 # lpm_extend_particle_array removed 412 415 # … … 715 718 716 719 CC = cc 717 CFLAGS = O 720 CFLAGS = O 718 721 719 722 F90 = … … 1628 1631 mod_kinds.o \ 1629 1632 modules.o \ 1633 netcdf_data_input_mod.o \ 1634 netcdf_interface_mod.o \ 1630 1635 surface_mod.o 1631 1636 swap_timelevel.o: \ 
palm/trunk/SOURCE/netcdf_interface_mod.f90
r3701 r3727 25 25 !  26 26 ! $Id$ 27 ! make several routines publicly available 28 ! 29 ! 3701 20190126 18:57:21Z knoop 27 30 ! Statement added to prevent compiler warning about unused variable 28 31 ! … … 599 602 END INTERFACE netcdf_create_file 600 603 604 INTERFACE netcdf_create_global_atts 605 MODULE PROCEDURE netcdf_create_global_atts 606 END INTERFACE netcdf_create_global_atts 607 601 608 INTERFACE netcdf_create_var 602 609 MODULE PROCEDURE netcdf_create_var … … 619 626 END INTERFACE netcdf_open_write_file 620 627 621 PUBLIC netcdf_create_file, netcdf_define_header, & 628 PUBLIC netcdf_create_att, netcdf_create_dim, netcdf_create_file, & 629 netcdf_create_global_atts, netcdf_create_var, netcdf_define_header, & 622 630 netcdf_handle_error, netcdf_open_write_file 623 631 … … 6756 6764 CALL message( 'define_netcdf_header', 'PA0259', 0, 0, 0, 6, 0 ) 6757 6765 6758 6766 6759 6767 CASE DEFAULT 6760 6768 
palm/trunk/SOURCE/surface_data_output_mod.f90
r3691 r3727 25 25 !  26 26 ! $Id$ 27 ! Enable NetCDF output for surface data (suehring, gronemeier) 28 ! 29 ! 3691 20190123 09:57:04Z suehring 27 30 ! Add output of surfaceparallel flow speed 28 31 ! … … 77 80 78 81 USE control_parameters, & 79 ONLY: surface_output 82 ONLY: coupling_char, data_output_during_spinup, end_time, & 83 message_string, run_description_header, simulated_time_at_begin, & 84 spinup_time, surface_output 80 85 81 86 USE grid_variables, & … … 85 90 ONLY: nxl, nxr, nys, nyn, nzb, nzt 86 91 92 #if defined( __netcdf ) 93 USE NETCDF 94 #endif 95 96 USE netcdf_data_input_mod, & 97 ONLY: init_model 98 99 USE netcdf_interface, & 100 ONLY: netcdf_create_att, netcdf_create_dim, netcdf_create_file, & 101 netcdf_create_global_atts, netcdf_create_var, netcdf_handle_error 102 87 103 USE pegrid 88 104 … … 100 116 INTEGER(iwp) :: npoints_total !< total number of points / vertices which define a surface element 101 117 118 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: s !< coordinate for NetCDF output, number of the surface element 119 102 120 REAL(wp) :: fillvalue = 9999.9_wp !< fillvalue for surface elements which are not defined 103 121 122 REAL(wp), DIMENSION(:), ALLOCATABLE :: azimuth !< azimuth orientation coordinate for NetCDF output 123 REAL(wp), DIMENSION(:), ALLOCATABLE :: es_utm !< EUTM coordinate for NetCDF output 124 REAL(wp), DIMENSION(:), ALLOCATABLE :: ns_utm !< EUTM coordinate for NetCDF output 125 REAL(wp), DIMENSION(:), ALLOCATABLE :: xs !< xcoordinate for NetCDF output 126 REAL(wp), DIMENSION(:), ALLOCATABLE :: ys !< ycoordinate for NetCDF output 127 REAL(wp), DIMENSION(:), ALLOCATABLE :: zs !< zcoordinate for NetCDF output 128 REAL(wp), DIMENSION(:), ALLOCATABLE :: zenith !< zenith orientation coordinate for NetCDF output 104 129 REAL(wp), DIMENSION(:), ALLOCATABLE :: var_out !< output variables 105 130 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: var_av !< variables used for averaging … … 109 134 110 135 CHARACTER(LEN=100), DIMENSION(300) :: data_output_surf = ' ' !< namelist variable which describes the output variables 111 CHARACTER(LEN=100), DIMENSION(0:1,300) :: dosurf = ' ' !< internal variable which describes the output variables and separates 112 !< averaged from nonaveraged output 136 CHARACTER(LEN=100), DIMENSION(0:1,300) :: dosurf = ' ' !< internal variable which describes the output variables 137 !< and separates averaged from nonaveraged output 138 CHARACTER(LEN=100), DIMENSION(0:1,300) :: dosurf_unit = ' ' !< internal variable which holds the unit of the given output variable 113 139 114 INTEGER(iwp) :: average_count_surf = 0 !< number of ensemble members used for averaging 115 INTEGER(iwp) :: dosurf_no(0:1) = 0 !< number of surface output quantities 140 INTEGER(iwp) :: average_count_surf = 0 !< number of ensemble members used for averaging 141 INTEGER(iwp) :: dosurf_no(0:1) = 0 !< number of surface output quantities 142 143 INTEGER(iwp) :: nc_stat !< error code for netcdf routines 144 INTEGER(iwp), SAVE :: oldmode !< save old setfillmode of netcdf file (not needed, but required for routine call) 145 146 INTEGER(iwp), DIMENSION(0:1) :: id_dim_s_surf !< netcdf ID for dimension s 147 INTEGER(iwp), DIMENSION(0:1) :: id_dim_time_surf !< netcdf ID for dimension time 148 INTEGER(iwp), DIMENSION(0:1) :: id_set_surf !< netcdf ID for file 149 INTEGER(iwp), DIMENSION(0:1) :: id_var_etum_surf !< netcdf ID for variable Es_UTM 150 INTEGER(iwp), DIMENSION(0:1) :: id_var_nutm_surf !< netcdf ID for variable Ns_UTM 151 INTEGER(iwp), DIMENSION(0:1) :: id_var_time_surf !< netcdf ID for variable time 152 INTEGER(iwp), DIMENSION(0:1) :: id_var_s_surf !< netcdf ID for variable s 153 INTEGER(iwp), DIMENSION(0:1) :: id_var_xs_surf !< netcdf ID for variable xs 154 INTEGER(iwp), DIMENSION(0:1) :: id_var_ys_surf !< netcdf ID for variable ys 155 INTEGER(iwp), DIMENSION(0:1) :: id_var_zs_surf !< netcdf ID for variable zs 156 INTEGER(iwp), DIMENSION(0:1) :: dosurf_time_count = 0 !< count of output time steps 157 INTEGER(iwp), DIMENSION(0:1) :: ntdim_surf !< number of output time steps 158 159 INTEGER(iwp), DIMENSION(0:1,300) :: id_var_dosurf !< netcdf ID for output variables 116 160 117 161 LOGICAL :: first_output(0:1) = .FALSE. !< true if first output was already called 162 LOGICAL :: to_netcdf = .TRUE. !< flag indicating parallel NetCDF output 163 LOGICAL :: to_vtk = .TRUE. !< flag indicating binary surfacedata output that can be further 164 !< processed to VTK format 118 165 119 166 REAL(wp) :: averaging_interval_surf = 9999999.9_wp !< averaging interval … … 125 172 REAL(wp) :: time_dosurf_av = 0.0_wp !< internal counter variable to check for averaged data output 126 173 127 128 174 TYPE(surf_out) :: surfaces !< variable which contains all required output information 129 175 … … 178 224 IMPLICIT NONE 179 225 226 CHARACTER (LEN=100) :: filename !< name of output file 227 CHARACTER (LEN=80) :: time_average_text !< string written to file attribute time_avg 228 CHARACTER (LEN=4000) :: var_list !< list of variables written to NetCDF file 229 230 INTEGER(iwp) :: av !< flag for averaged (=1) and nonaveraged (=0) data 180 231 INTEGER(iwp) :: i !< grid index in xdirection, also running variable for counting nonaverage data output 181 232 INTEGER(iwp) :: j !< grid index in ydirection, also running variable for counting average data output … … 183 234 INTEGER(iwp) :: l !< running index for surfaceelement orientation 184 235 INTEGER(iwp) :: m !< running index for surface elements 236 INTEGER(iwp) :: mm !< local counting variable for surface elements 185 237 INTEGER(iwp) :: npg !< counter variable for all surface elements ( or polygons ) 186 238 INTEGER(iwp) :: point_index_count !< local counter variable for point index 239 INTEGER(iwp) :: start_count !< local start counter for the surface index 187 240 188 241 INTEGER(iwp), DIMENSION(0:numprocs1) :: num_points_on_pe !< array which contains the number of points on all mpi ranks 242 INTEGER(iwp), DIMENSION(0:numprocs1) :: num_surfaces_on_pe !< array which contains the number of surfaces on all mpi ranks 189 243 INTEGER(iwp), ALLOCATABLE, DIMENSION(:,:,:) :: point_index !< dummy array used to check where the reference points for surface polygons are located 190 244 245 REAL(wp) :: az !< azimuth angle, indicated the vertical orientation of a surface element 246 REAL(wp) :: off_x !< grid offset in xdirection between the stored grid index and the actual wall 247 REAL(wp) :: off_y !< grid offset in ydirection between the stored grid index and the actual wall 248 249 REAL(wp), DIMENSION(:), ALLOCATABLE :: netcdf_data_1d !< dummy array to output 1D data into netcdf file 250 251 191 252 ! 192 253 ! Determine the number of surface elements on subdomain … … 198 259 + surf_def_v(3)%ns + surf_lsm_v(3)%ns + surf_usm_v(3)%ns !eastwardfacing 199 260 ! 261 ! Determine the total number of surfaces in the model domain 262 #if defined( __parallel ) 263 CALL MPI_ALLREDUCE( surfaces%ns, surfaces%ns_total, 1, & 264 MPI_INTEGER, MPI_SUM, comm2d, ierr ) 265 #else 266 surfaces%ns_total = surfaces%ns 267 #endif 268 ! 200 269 ! Allocate output variable and set to _FillValue attribute 201 270 ALLOCATE ( surfaces%var_out(1:surfaces%ns) ) … … 209 278 ENDIF 210 279 ! 211 ! I nitialize pointsand polygon data.280 ! If output to VTK format is enabled, initialize point and polygon data. 212 281 ! In a first step, count the number of points which are defining 213 282 ! the surfaces and the polygons. 214 ALLOCATE( point_index(nzb:nzt+1,nys1:nyn+1,nxl1:nxr+1) ) 215 point_index = 1 216 ! 217 ! Horizontal default surfaces 218 surfaces%npoints = 0 219 DO l = 0, 1 220 DO m = 1, surf_def_h(0)%ns 221 ! 222 ! Determine the indices of the respective grid cell inside the topography 223 i = surf_def_h(0)%i(m) + surf_def_h(0)%ioff 224 j = surf_def_h(0)%j(m) + surf_def_h(0)%joff 225 k = surf_def_h(0)%k(m) + surf_def_h(0)%koff 226 ! 227 ! Check if the vertices that define the surface element are already 228 ! defined, if not, increment the counter. 283 IF ( to_vtk ) THEN 284 ALLOCATE( point_index(nzb:nzt+1,nys1:nyn+1,nxl1:nxr+1) ) 285 point_index = 1 286 ! 287 ! Horizontal default surfaces 288 surfaces%npoints = 0 289 DO l = 0, 1 290 DO m = 1, surf_def_h(0)%ns 291 ! 292 ! Determine the indices of the respective grid cell inside the topography 293 i = surf_def_h(0)%i(m) + surf_def_h(0)%ioff 294 j = surf_def_h(0)%j(m) + surf_def_h(0)%joff 295 k = surf_def_h(0)%k(m) + surf_def_h(0)%koff 296 ! 297 ! Check if the vertices that define the surface element are already 298 ! defined, if not, increment the counter. 299 IF ( point_index(k,j,i) < 0 ) THEN 300 surfaces%npoints = surfaces%npoints + 1 301 point_index(k,j,i) = surfaces%npoints  1 302 ENDIF 303 IF ( point_index(k,j,i+1) < 0 ) THEN 304 surfaces%npoints = surfaces%npoints + 1 305 point_index(k,j,i+1) = surfaces%npoints  1 306 ENDIF 307 IF ( point_index(k,j+1,i+1) < 0 ) THEN 308 surfaces%npoints = surfaces%npoints + 1 309 point_index(k,j+1,i+1) = surfaces%npoints  1 310 ENDIF 311 IF ( point_index(k,j+1,i) < 0 ) THEN 312 surfaces%npoints = surfaces%npoints + 1 313 point_index(k,j+1,i) = surfaces%npoints  1 314 ENDIF 315 ENDDO 316 ENDDO 317 DO m = 1, surf_lsm_h%ns 318 i = surf_lsm_h%i(m) + surf_lsm_h%ioff 319 j = surf_lsm_h%j(m) + surf_lsm_h%joff 320 k = surf_lsm_h%k(m) + surf_lsm_h%koff 321 229 322 IF ( point_index(k,j,i) < 0 ) THEN 230 323 surfaces%npoints = surfaces%npoints + 1 … … 242 335 surfaces%npoints = surfaces%npoints + 1 243 336 point_index(k,j+1,i) = surfaces%npoints  1 244 ENDIF 337 ENDIF 245 338 ENDDO 246 ENDDO 247 DO m = 1, surf_lsm_h%ns 248 i = surf_lsm_h%i(m) + surf_lsm_h%ioff 249 j = surf_lsm_h%j(m) + surf_lsm_h%joff 250 k = surf_lsm_h%k(m) + surf_lsm_h%koff 251 252 IF ( point_index(k,j,i) < 0 ) THEN 253 surfaces%npoints = surfaces%npoints + 1 254 point_index(k,j,i) = surfaces%npoints  1 255 ENDIF 256 IF ( point_index(k,j,i+1) < 0 ) THEN 257 surfaces%npoints = surfaces%npoints + 1 258 point_index(k,j,i+1) = surfaces%npoints  1 259 ENDIF 260 IF ( point_index(k,j+1,i+1) < 0 ) THEN 261 surfaces%npoints = surfaces%npoints + 1 262 point_index(k,j+1,i+1) = surfaces%npoints  1 263 ENDIF 264 IF ( point_index(k,j+1,i) < 0 ) THEN 265 surfaces%npoints = surfaces%npoints + 1 266 point_index(k,j+1,i) = surfaces%npoints  1 267 ENDIF 268 ENDDO 269 DO m = 1, surf_usm_h%ns 270 i = surf_usm_h%i(m) + surf_usm_h%ioff 271 j = surf_usm_h%j(m) + surf_usm_h%joff 272 k = surf_usm_h%k(m) + surf_usm_h%koff 273 274 IF ( point_index(k,j,i) < 0 ) THEN 275 surfaces%npoints = surfaces%npoints + 1 276 point_index(k,j,i) = surfaces%npoints  1 277 ENDIF 278 IF ( point_index(k,j,i+1) < 0 ) THEN 279 surfaces%npoints = surfaces%npoints + 1 280 point_index(k,j,i+1) = surfaces%npoints  1 281 ENDIF 282 IF ( point_index(k,j+1,i+1) < 0 ) THEN 283 surfaces%npoints = surfaces%npoints + 1 284 point_index(k,j+1,i+1) = surfaces%npoints  1 285 ENDIF 286 IF ( point_index(k,j+1,i) < 0 ) THEN 287 surfaces%npoints = surfaces%npoints + 1 288 point_index(k,j+1,i) = surfaces%npoints  1 289 ENDIF 290 ENDDO 291 ! 292 ! Vertical surfaces 293 DO l = 0, 3 294 DO m = 1, surf_def_v(l)%ns 295 ! 296 ! Determine the indices of the respective grid cell inside the topography. 297 ! Please note, jindex for northfacing surfaces ( l==0 ) is 298 ! identical to the reference jindex outside the grid box. 299 ! Equivalent for eastfacing surfaces and iindex. 300 i = surf_def_v(l)%i(m) + MERGE( surf_def_v(l)%ioff, 0, l == 3 ) 301 j = surf_def_v(l)%j(m) + MERGE( surf_def_v(l)%joff, 0, l == 1 ) 302 k = surf_def_v(l)%k(m) + surf_def_v(l)%koff 303 ! 304 ! Lower left /front vertex 339 DO m = 1, surf_usm_h%ns 340 i = surf_usm_h%i(m) + surf_usm_h%ioff 341 j = surf_usm_h%j(m) + surf_usm_h%joff 342 k = surf_usm_h%k(m) + surf_usm_h%koff 343 305 344 IF ( point_index(k,j,i) < 0 ) THEN 306 345 surfaces%npoints = surfaces%npoints + 1 307 346 point_index(k,j,i) = surfaces%npoints  1 308 ENDIF 309 ! 310 ! Upper / lower right index for north and southfacing surfaces 311 IF ( l == 0 .OR. l == 1 ) THEN 347 ENDIF 348 IF ( point_index(k,j,i+1) < 0 ) THEN 349 surfaces%npoints = surfaces%npoints + 1 350 point_index(k,j,i+1) = surfaces%npoints  1 351 ENDIF 352 IF ( point_index(k,j+1,i+1) < 0 ) THEN 353 surfaces%npoints = surfaces%npoints + 1 354 point_index(k,j+1,i+1) = surfaces%npoints  1 355 ENDIF 356 IF ( point_index(k,j+1,i) < 0 ) THEN 357 surfaces%npoints = surfaces%npoints + 1 358 point_index(k,j+1,i) = surfaces%npoints  1 359 ENDIF 360 ENDDO 361 ! 362 ! Vertical surfaces 363 DO l = 0, 3 364 DO m = 1, surf_def_v(l)%ns 365 ! 366 ! Determine the indices of the respective grid cell inside the 367 ! topography. Please note, jindex for northfacing surfaces 368 ! ( l==0 ) is identical to the reference jindex outside the grid 369 ! box. Equivalent for eastfacing surfaces and iindex. 370 i = surf_def_v(l)%i(m) + MERGE( surf_def_v(l)%ioff, 0, l == 3 ) 371 j = surf_def_v(l)%j(m) + MERGE( surf_def_v(l)%joff, 0, l == 1 ) 372 k = surf_def_v(l)%k(m) + surf_def_v(l)%koff 373 ! 374 ! Lower left /front vertex 375 IF ( point_index(k,j,i) < 0 ) THEN 376 surfaces%npoints = surfaces%npoints + 1 377 point_index(k,j,i) = surfaces%npoints  1 378 ENDIF 379 ! 380 ! Upper / lower right index for north and southfacing surfaces 381 IF ( l == 0 .OR. l == 1 ) THEN 382 IF ( point_index(k,j,i+1) < 0 ) THEN 383 surfaces%npoints = surfaces%npoints + 1 384 point_index(k,j,i+1) = surfaces%npoints  1 385 ENDIF 386 IF ( point_index(k+1,j,i+1) < 0 ) THEN 387 surfaces%npoints = surfaces%npoints + 1 388 point_index(k+1,j,i+1) = surfaces%npoints  1 389 ENDIF 390 ! 391 ! Upper / lower front index for east and westfacing surfaces 392 ELSEIF ( l == 2 .OR. l == 3 ) THEN 393 IF ( point_index(k,j+1,i) < 0 ) THEN 394 surfaces%npoints = surfaces%npoints + 1 395 point_index(k,j+1,i) = surfaces%npoints  1 396 ENDIF 397 IF ( point_index(k+1,j+1,i) < 0 ) THEN 398 surfaces%npoints = surfaces%npoints + 1 399 point_index(k+1,j+1,i) = surfaces%npoints  1 400 ENDIF 401 ENDIF 402 ! 403 ! Upper left / front vertex 404 IF ( point_index(k+1,j,i) < 0 ) THEN 405 surfaces%npoints = surfaces%npoints + 1 406 point_index(k+1,j,i) = surfaces%npoints  1 407 ENDIF 408 ENDDO 409 DO m = 1, surf_lsm_v(l)%ns 410 i = surf_lsm_v(l)%i(m) + MERGE( surf_lsm_v(l)%ioff, 0, l == 3 ) 411 j = surf_lsm_v(l)%j(m) + MERGE( surf_lsm_v(l)%joff, 0, l == 1 ) 412 k = surf_lsm_v(l)%k(m) + surf_lsm_v(l)%koff 413 ! 414 ! Lower left /front vertex 415 IF ( point_index(k,j,i) < 0 ) THEN 416 surfaces%npoints = surfaces%npoints + 1 417 point_index(k,j,i) = surfaces%npoints  1 418 ENDIF 419 ! 420 ! Upper / lower right index for north and southfacing surfaces 421 IF ( l == 0 .OR. l == 1 ) THEN 422 IF ( point_index(k,j,i+1) < 0 ) THEN 423 surfaces%npoints = surfaces%npoints + 1 424 point_index(k,j,i+1) = surfaces%npoints  1 425 ENDIF 426 IF ( point_index(k+1,j,i+1) < 0 ) THEN 427 surfaces%npoints = surfaces%npoints + 1 428 point_index(k+1,j,i+1) = surfaces%npoints  1 429 ENDIF 430 ! 431 ! Upper / lower front index for east and westfacing surfaces 432 ELSEIF ( l == 2 .OR. l == 3 ) THEN 433 IF ( point_index(k,j+1,i) < 0 ) THEN 434 surfaces%npoints = surfaces%npoints + 1 435 point_index(k,j+1,i) = surfaces%npoints  1 436 ENDIF 437 IF ( point_index(k+1,j+1,i) < 0 ) THEN 438 surfaces%npoints = surfaces%npoints + 1 439 point_index(k+1,j+1,i) = surfaces%npoints  1 440 ENDIF 441 ENDIF 442 ! 443 ! Upper left / front vertex 444 IF ( point_index(k+1,j,i) < 0 ) THEN 445 surfaces%npoints = surfaces%npoints + 1 446 point_index(k+1,j,i) = surfaces%npoints  1 447 ENDIF 448 ENDDO 449 450 DO m = 1, surf_usm_v(l)%ns 451 i = surf_usm_v(l)%i(m) + MERGE( surf_usm_v(l)%ioff, 0, l == 3 ) 452 j = surf_usm_v(l)%j(m) + MERGE( surf_usm_v(l)%joff, 0, l == 1 ) 453 k = surf_usm_v(l)%k(m) + surf_usm_v(l)%koff 454 ! 455 ! Lower left /front vertex 456 IF ( point_index(k,j,i) < 0 ) THEN 457 surfaces%npoints = surfaces%npoints + 1 458 point_index(k,j,i) = surfaces%npoints  1 459 ENDIF 460 ! 461 ! Upper / lower right index for north and southfacing surfaces 462 IF ( l == 0 .OR. l == 1 ) THEN 463 IF ( point_index(k,j,i+1) < 0 ) THEN 464 surfaces%npoints = surfaces%npoints + 1 465 point_index(k,j,i+1) = surfaces%npoints  1 466 ENDIF 467 IF ( point_index(k+1,j,i+1) < 0 ) THEN 468 surfaces%npoints = surfaces%npoints + 1 469 point_index(k+1,j,i+1) = surfaces%npoints  1 470 ENDIF 471 ! 472 ! Upper / lower front index for east and westfacing surfaces 473 ELSEIF ( l == 2 .OR. l == 3 ) THEN 474 IF ( point_index(k,j+1,i) < 0 ) THEN 475 surfaces%npoints = surfaces%npoints + 1 476 point_index(k,j+1,i) = surfaces%npoints  1 477 ENDIF 478 IF ( point_index(k+1,j+1,i) < 0 ) THEN 479 surfaces%npoints = surfaces%npoints + 1 480 point_index(k+1,j+1,i) = surfaces%npoints  1 481 ENDIF 482 ENDIF 483 ! 484 ! Upper left / front vertex 485 IF ( point_index(k+1,j,i) < 0 ) THEN 486 surfaces%npoints = surfaces%npoints + 1 487 point_index(k+1,j,i) = surfaces%npoints  1 488 ENDIF 489 ENDDO 490 491 ENDDO 492 ! 493 ! Allocate the number of points and polygons. Note, the number of polygons 494 ! is identical to the number of surfaces elements, whereas the number 495 ! of points (vertices), which define the polygons, can be larger. 496 ALLOCATE( surfaces%points(3,1:surfaces%npoints) ) 497 ALLOCATE( surfaces%polygons(5,1:surfaces%ns) ) 498 ! 499 ! Note, PARAVIEW expects consecutively ordered points, in order to 500 ! unambiguously identify surfaces. 501 ! Hence, all PEs should know where they start counting, depending on the 502 ! number of points on the other PE's with lower MPI rank. 503 #if defined( __parallel ) 504 CALL MPI_ALLGATHER( surfaces%npoints, 1, MPI_INTEGER, & 505 num_points_on_pe, 1, MPI_INTEGER, comm2d, ierr ) 506 #else 507 num_points_on_pe = surfaces%npoints 508 #endif 509 510 ! 511 ! After the number of vertices is counted, repeat the loops and define the 512 ! vertices. Start with the horizontal default surfaces. 513 ! First, however, determine the offset where couting of points should be 514 ! started, which is the sum of points of all PE's with lower MPI rank. 515 i = 0 516 point_index_count = 0 517 DO WHILE ( i < myid .AND. i <= SIZE( num_points_on_pe ) ) 518 point_index_count = point_index_count + num_points_on_pe(i) 519 i = i + 1 520 ENDDO 521 522 surfaces%npoints = 0 523 point_index = 1 524 npg = 0 525 526 DO l = 0, 1 527 DO m = 1, surf_def_h(0)%ns 528 ! 529 ! Determine the indices of the respective grid cell inside the 530 ! topography. 531 i = surf_def_h(0)%i(m) + surf_def_h(0)%ioff 532 j = surf_def_h(0)%j(m) + surf_def_h(0)%joff 533 k = surf_def_h(0)%k(m) + surf_def_h(0)%koff 534 ! 535 ! Check if the vertices that define the surface element are 536 ! already defined, if not, increment the counter. 537 IF ( point_index(k,j,i) < 0 ) THEN 538 surfaces%npoints = surfaces%npoints + 1 539 point_index(k,j,i) = point_index_count 540 point_index_count = point_index_count + 1 541 surfaces%points(1,surfaces%npoints) = ( i  0.5_wp ) * dx 542 surfaces%points(2,surfaces%npoints) = ( j  0.5_wp ) * dy 543 surfaces%points(3,surfaces%npoints) = zw(k) 544 ENDIF 312 545 IF ( point_index(k,j,i+1) < 0 ) THEN 313 546 surfaces%npoints = surfaces%npoints + 1 314 point_index(k,j,i+1) = surfaces%npoints  1 315 ENDIF 316 IF ( point_index(k+1,j,i+1) < 0 ) THEN 547 point_index(k,j,i+1) = point_index_count 548 point_index_count = point_index_count + 1 549 surfaces%points(1,surfaces%npoints) = ( i + 1  0.5_wp ) * dx 550 surfaces%points(2,surfaces%npoints) = ( j  0.5_wp ) * dy 551 surfaces%points(3,surfaces%npoints) = zw(k) 552 ENDIF 553 IF ( point_index(k,j+1,i+1) < 0 ) THEN 317 554 surfaces%npoints = surfaces%npoints + 1 318 point_index(k+1,j,i+1) = surfaces%npoints  1 319 ENDIF 320 ! 321 ! Upper / lower front index for east and westfacing surfaces 322 ELSEIF ( l == 2 .OR. l == 3 ) THEN 555 point_index(k,j+1,i+1) = point_index_count 556 point_index_count = point_index_count + 1 557 surfaces%points(1,surfaces%npoints) = ( i + 1  0.5_wp ) * dx 558 surfaces%points(2,surfaces%npoints) = ( j + 1  0.5_wp ) * dy 559 surfaces%points(3,surfaces%npoints) = zw(k) 560 ENDIF 323 561 IF ( point_index(k,j+1,i) < 0 ) THEN 324 562 surfaces%npoints = surfaces%npoints + 1 325 point_index(k,j+1,i) = surfaces%npoints  1 326 ENDIF 327 IF ( point_index(k+1,j+1,i) < 0 ) THEN 328 surfaces%npoints = surfaces%npoints + 1 329 point_index(k+1,j+1,i) = surfaces%npoints  1 330 ENDIF 331 ENDIF 332 ! 333 ! Upper left / front vertex 334 IF ( point_index(k+1,j,i) < 0 ) THEN 335 surfaces%npoints = surfaces%npoints + 1 336 point_index(k+1,j,i) = surfaces%npoints  1 337 ENDIF 563 point_index(k,j+1,i) = point_index_count 564 point_index_count = point_index_count + 1 565 surfaces%points(1,surfaces%npoints) = ( i  0.5_wp ) * dx 566 surfaces%points(2,surfaces%npoints) = ( j + 1  0.5_wp ) * dy 567 surfaces%points(3,surfaces%npoints) = zw(k) 568 ENDIF 569 570 npg = npg + 1 571 surfaces%polygons(1,npg) = 4 572 surfaces%polygons(2,npg) = point_index(k,j,i) 573 surfaces%polygons(3,npg) = point_index(k,j,i+1) 574 surfaces%polygons(4,npg) = point_index(k,j+1,i+1) 575 surfaces%polygons(5,npg) = point_index(k,j+1,i) 576 ENDDO 338 577 ENDDO 339 DO m = 1, surf_lsm_v(l)%ns 340 i = surf_lsm_v(l)%i(m) + MERGE( surf_lsm_v(l)%ioff, 0, l == 3 ) 341 j = surf_lsm_v(l)%j(m) + MERGE( surf_lsm_v(l)%joff, 0, l == 1 ) 342 k = surf_lsm_v(l)%k(m) + surf_lsm_v(l)%koff 343 ! 344 ! Lower left /front vertex 345 IF ( point_index(k,j,i) < 0 ) THEN 346 surfaces%npoints = surfaces%npoints + 1 347 point_index(k,j,i) = surfaces%npoints  1 348 ENDIF 349 ! 350 ! Upper / lower right index for north and southfacing surfaces 351 IF ( l == 0 .OR. l == 1 ) THEN 352 IF ( point_index(k,j,i+1) < 0 ) THEN 353 surfaces%npoints = surfaces%npoints + 1 354 point_index(k,j,i+1) = surfaces%npoints  1 355 ENDIF 356 IF ( point_index(k+1,j,i+1) < 0 ) THEN 357 surfaces%npoints = surfaces%npoints + 1 358 point_index(k+1,j,i+1) = surfaces%npoints  1 359 ENDIF 360 ! 361 ! Upper / lower front index for east and westfacing surfaces 362 ELSEIF ( l == 2 .OR. l == 3 ) THEN 363 IF ( point_index(k,j+1,i) < 0 ) THEN 364 surfaces%npoints = surfaces%npoints + 1 365 point_index(k,j+1,i) = surfaces%npoints  1 366 ENDIF 367 IF ( point_index(k+1,j+1,i) < 0 ) THEN 368 surfaces%npoints = surfaces%npoints + 1 369 point_index(k+1,j+1,i) = surfaces%npoints  1 370 ENDIF 371 ENDIF 372 ! 373 ! Upper left / front vertex 374 IF ( point_index(k+1,j,i) < 0 ) THEN 375 surfaces%npoints = surfaces%npoints + 1 376 point_index(k+1,j,i) = surfaces%npoints  1 377 ENDIF 378 ENDDO 379 380 DO m = 1, surf_usm_v(l)%ns 381 i = surf_usm_v(l)%i(m) + MERGE( surf_usm_v(l)%ioff, 0, l == 3 ) 382 j = surf_usm_v(l)%j(m) + MERGE( surf_usm_v(l)%joff, 0, l == 1 ) 383 k = surf_usm_v(l)%k(m) + surf_usm_v(l)%koff 384 ! 385 ! Lower left /front vertex 386 IF ( point_index(k,j,i) < 0 ) THEN 387 surfaces%npoints = surfaces%npoints + 1 388 point_index(k,j,i) = surfaces%npoints  1 389 ENDIF 390 ! 391 ! Upper / lower right index for north and southfacing surfaces 392 IF ( l == 0 .OR. l == 1 ) THEN 393 IF ( point_index(k,j,i+1) < 0 ) THEN 394 surfaces%npoints = surfaces%npoints + 1 395 point_index(k,j,i+1) = surfaces%npoints  1 396 ENDIF 397 IF ( point_index(k+1,j,i+1) < 0 ) THEN 398 surfaces%npoints = surfaces%npoints + 1 399 point_index(k+1,j,i+1) = surfaces%npoints  1 400 ENDIF 401 ! 402 ! Upper / lower front index for east and westfacing surfaces 403 ELSEIF ( l == 2 .OR. l == 3 ) THEN 404 IF ( point_index(k,j+1,i) < 0 ) THEN 405 surfaces%npoints = surfaces%npoints + 1 406 point_index(k,j+1,i) = surfaces%npoints  1 407 ENDIF 408 IF ( point_index(k+1,j+1,i) < 0 ) THEN 409 surfaces%npoints = surfaces%npoints + 1 410 point_index(k+1,j+1,i) = surfaces%npoints  1 411 ENDIF 412 ENDIF 413 ! 414 ! Upper left / front vertex 415 IF ( point_index(k+1,j,i) < 0 ) THEN 416 surfaces%npoints = surfaces%npoints + 1 417 point_index(k+1,j,i) = surfaces%npoints  1 418 ENDIF 419 ENDDO 420 421 ENDDO 422 ! 423 ! Allocate the number of points and polygons. Note, the number of polygons 424 ! is identical to the number of surfaces elements, whereas the number 425 ! of points (vertices) which define the polygons can be larger. 426 ALLOCATE( surfaces%points(3,1:surfaces%npoints) ) 427 ALLOCATE( surfaces%polygons(5,1:surfaces%ns) ) 428 ! 429 ! Note, PARAVIEW expects consecutively ordered points which can be 430 ! unambiguously identified. 431 ! Hence, all PEs know where they should start counting, depending on the 432 ! number of points on the other PE's with lower MPI rank. 433 #if defined( __parallel ) 434 CALL MPI_ALLGATHER( surfaces%npoints, 1, MPI_INTEGER, & 435 num_points_on_pe, 1, MPI_INTEGER, comm2d, ierr ) 436 #else 437 num_points_on_pe = surfaces%npoints 438 #endif 439 440 ! 441 ! After the number of vertices is counted, repeat the loops and define the 442 ! vertices. Start with the horizontal default surfaces 443 ! First, however, determine the offset where couting of points should be 444 ! started, which is the sum of points of all PE's with lower MPI rank. 445 i = 0 446 point_index_count = 0 447 DO WHILE ( i < myid .AND. i <= SIZE( num_points_on_pe ) ) 448 point_index_count = point_index_count + num_points_on_pe(i) 449 i = i + 1 450 ENDDO 451 452 453 surfaces%npoints = 0 454 point_index = 1 455 npg = 0 456 457 DO l = 0, 1 458 DO m = 1, surf_def_h(0)%ns 459 ! 460 ! Determine the indices of the respective grid cell inside the topography 461 i = surf_def_h(0)%i(m) + surf_def_h(0)%ioff 462 j = surf_def_h(0)%j(m) + surf_def_h(0)%joff 463 k = surf_def_h(0)%k(m) + surf_def_h(0)%koff 464 ! 465 ! Check if the vertices that define the surface element are already 466 ! defined, if not, increment the counter. 578 DO m = 1, surf_lsm_h%ns 579 i = surf_lsm_h%i(m) + surf_lsm_h%ioff 580 j = surf_lsm_h%j(m) + surf_lsm_h%joff 581 k = surf_lsm_h%k(m) + surf_lsm_h%koff 467 582 IF ( point_index(k,j,i) < 0 ) THEN 468 583 surfaces%npoints = surfaces%npoints + 1 469 584 point_index(k,j,i) = point_index_count 470 point_index_count = point_index_count + 1 585 point_index_count = point_index_count + 1 471 586 surfaces%points(1,surfaces%npoints) = ( i  0.5_wp ) * dx 472 587 surfaces%points(2,surfaces%npoints) = ( j  0.5_wp ) * dy … … 503 618 surfaces%polygons(3,npg) = point_index(k,j,i+1) 504 619 surfaces%polygons(4,npg) = point_index(k,j+1,i+1) 505 surfaces%polygons(5,npg) = point_index(k,j+1,i) 620 surfaces%polygons(5,npg) = point_index(k,j+1,i) 506 621 ENDDO 507 ENDDO 508 DO m = 1, surf_lsm_h%ns 509 i = surf_lsm_h%i(m) + surf_lsm_h%ioff 510 j = surf_lsm_h%j(m) + surf_lsm_h%joff 511 k = surf_lsm_h%k(m) + surf_lsm_h%koff 512 IF ( point_index(k,j,i) < 0 ) THEN 513 surfaces%npoints = surfaces%npoints + 1 514 point_index(k,j,i) = point_index_count 515 point_index_count = point_index_count + 1 516 surfaces%points(1,surfaces%npoints) = ( i  0.5_wp ) * dx 517 surfaces%points(2,surfaces%npoints) = ( j  0.5_wp ) * dy 518 surfaces%points(3,surfaces%npoints) = zw(k) 519 ENDIF 520 IF ( point_index(k,j,i+1) < 0 ) THEN 521 surfaces%npoints = surfaces%npoints + 1 522 point_index(k,j,i+1) = point_index_count 523 point_index_count = point_index_count + 1 524 surfaces%points(1,surfaces%npoints) = ( i + 1  0.5_wp ) * dx 525 surfaces%points(2,surfaces%npoints) = ( j  0.5_wp ) * dy 526 surfaces%points(3,surfaces%npoints) = zw(k) 527 ENDIF 528 IF ( point_index(k,j+1,i+1) < 0 ) THEN 529 surfaces%npoints = surfaces%npoints + 1 530 point_index(k,j+1,i+1) = point_index_count 531 point_index_count = point_index_count + 1 532 surfaces%points(1,surfaces%npoints) = ( i + 1  0.5_wp ) * dx 533 surfaces%points(2,surfaces%npoints) = ( j + 1  0.5_wp ) * dy 534 surfaces%points(3,surfaces%npoints) = zw(k) 535 ENDIF 536 IF ( point_index(k,j+1,i) < 0 ) THEN 537 surfaces%npoints = surfaces%npoints + 1 538 point_index(k,j+1,i) = point_index_count 539 point_index_count = point_index_count + 1 540 surfaces%points(1,surfaces%npoints) = ( i  0.5_wp ) * dx 541 surfaces%points(2,surfaces%npoints) = ( j + 1  0.5_wp ) * dy 542 surfaces%points(3,surfaces%npoints) = zw(k) 543 ENDIF 544 545 npg = npg + 1 546 surfaces%polygons(1,npg) = 4 547 surfaces%polygons(2,npg) = point_index(k,j,i) 548 surfaces%polygons(3,npg) = point_index(k,j,i+1) 549 surfaces%polygons(4,npg) = point_index(k,j+1,i+1) 550 surfaces%polygons(5,npg) = point_index(k,j+1,i) 551 ENDDO 552 553 DO m = 1, surf_usm_h%ns 554 i = surf_usm_h%i(m) + surf_usm_h%ioff 555 j = surf_usm_h%j(m) + surf_usm_h%joff 556 k = surf_usm_h%k(m) + surf_usm_h%koff 557 558 IF ( point_index(k,j,i) < 0 ) THEN 559 surfaces%npoints = surfaces%npoints + 1 560 point_index(k,j,i) = point_index_count 561 point_index_count = point_index_count + 1 562 surfaces%points(1,surfaces%npoints) = ( i  0.5_wp ) * dx 563 surfaces%points(2,surfaces%npoints) = ( j  0.5_wp ) * dy 564 surfaces%points(3,surfaces%npoints) = zw(k) 565 ENDIF 566 IF ( point_index(k,j,i+1) < 0 ) THEN 567 surfaces%npoints = surfaces%npoints + 1 568 point_index(k,j,i+1) = point_index_count 569 point_index_count = point_index_count + 1 570 surfaces%points(1,surfaces%npoints) = ( i + 1  0.5_wp ) * dx 571 surfaces%points(2,surfaces%npoints) = ( j  0.5_wp ) * dy 572 surfaces%points(3,surfaces%npoints) = zw(k) 573 ENDIF 574 IF ( point_index(k,j+1,i+1) < 0 ) THEN 575 surfaces%npoints = surfaces%npoints + 1 576 point_index(k,j+1,i+1) = point_index_count 577 point_index_count = point_index_count + 1 578 surfaces%points(1,surfaces%npoints) = ( i + 1  0.5_wp ) * dx 579 surfaces%points(2,surfaces%npoints) = ( j + 1  0.5_wp ) * dy 580 surfaces%points(3,surfaces%npoints) = zw(k) 581 ENDIF 582 IF ( point_index(k,j+1,i) < 0 ) THEN 583 surfaces%npoints = surfaces%npoints + 1 584 point_index(k,j+1,i) = point_index_count 585 point_index_count = point_index_count + 1 586 surfaces%points(1,surfaces%npoints) = ( i  0.5_wp ) * dx 587 surfaces%points(2,surfaces%npoints) = ( j + 1  0.5_wp ) * dy 588 surfaces%points(3,surfaces%npoints) = zw(k) 589 ENDIF 590 591 npg = npg + 1 592 surfaces%polygons(1,npg) = 4 593 surfaces%polygons(2,npg) = point_index(k,j,i) 594 surfaces%polygons(3,npg) = point_index(k,j,i+1) 595 surfaces%polygons(4,npg) = point_index(k,j+1,i+1) 596 surfaces%polygons(5,npg) = point_index(k,j+1,i) 597 ENDDO 598 599 DO l = 0, 3 600 DO m = 1, surf_def_v(l)%ns 601 ! 602 ! Determine the indices of the respective grid cell inside the topography. 603 ! Please note, jindex for northfacing surfaces ( l==0 ) is 604 ! identical to the reference jindex outside the grid box. 605 ! Equivalent for eastfacing surfaces and iindex. 606 i = surf_def_v(l)%i(m) + MERGE( surf_def_v(l)%ioff, 0, l == 3 ) 607 j = surf_def_v(l)%j(m) + MERGE( surf_def_v(l)%joff, 0, l == 1 ) 608 k = surf_def_v(l)%k(m) + surf_def_v(l)%koff 609 ! 610 ! Lower left /front vertex 622 623 DO m = 1, surf_usm_h%ns 624 i = surf_usm_h%i(m) + surf_usm_h%ioff 625 j = surf_usm_h%j(m) + surf_usm_h%joff 626 k = surf_usm_h%k(m) + surf_usm_h%koff 627 611 628 IF ( point_index(k,j,i) < 0 ) THEN 612 629 surfaces%npoints = surfaces%npoints + 1 … … 615 632 surfaces%points(1,surfaces%npoints) = ( i  0.5_wp ) * dx 616 633 surfaces%points(2,surfaces%npoints) = ( j  0.5_wp ) * dy 617 surfaces%points(3,surfaces%npoints) = zw(k1) 618 ENDIF 619 ! 620 ! Upper / lower right index for north and southfacing surfaces 621 IF ( l == 0 .OR. l == 1 ) THEN 622 IF ( point_index(k,j,i+1) < 0 ) THEN 623 surfaces%npoints = surfaces%npoints + 1 624 point_index(k,j,i+1) = point_index_count 625 point_index_count = point_index_count + 1 626 surfaces%points(1,surfaces%npoints) = ( i + 1  0.5_wp ) * dx 627 surfaces%points(2,surfaces%npoints) = ( j  0.5_wp ) * dy 628 surfaces%points(3,surfaces%npoints) = zw(k1) 629 ENDIF 630 IF ( point_index(k+1,j,i+1) < 0 ) THEN 631 surfaces%npoints = surfaces%npoints + 1 632 point_index(k+1,j,i+1) = point_index_count 633 point_index_count = point_index_count + 1 634 surfaces%points(1,surfaces%npoints) = ( i + 1  0.5_wp ) * dx 635 surfaces%points(2,surfaces%npoints) = ( j  0.5_wp ) * dy 636 surfaces%points(3,surfaces%npoints) = zw(k) 637 ENDIF 638 ! 639 ! Upper / lower front index for east and westfacing surfaces 640 ELSEIF ( l == 2 .OR. l == 3 ) THEN 641 IF ( point_index(k,j+1,i) < 0 ) THEN 642 surfaces%npoints = surfaces%npoints + 1 643 point_index(k,j+1,i) = point_index_count 644 point_index_count = point_index_count + 1 645 surfaces%points(1,surfaces%npoints) = ( i  0.5_wp ) * dx 646 surfaces%points(2,surfaces%npoints) = ( j + 1  0.5_wp ) * dy 647 surfaces%points(3,surfaces%npoints) = zw(k1) 648 ENDIF 649 IF ( point_index(k+1,j+1,i) < 0 ) THEN 650 surfaces%npoints = surfaces%npoints + 1 651 point_index(k+1,j+1,i) = point_index_count 652 point_index_count = point_index_count + 1 653 surfaces%points(1,surfaces%npoints) = ( i  0.5_wp ) * dx 654 surfaces%points(2,surfaces%npoints) = ( j + 1  0.5_wp ) * dy 655 surfaces%points(3,surfaces%npoints) = zw(k) 656 ENDIF 634 surfaces%points(3,surfaces%npoints) = zw(k) 657 635 ENDIF 658 ! 659 ! Upper left / front vertex 660 IF ( point_index(k+1,j,i) < 0 ) THEN 636 IF ( point_index(k,j,i+1) < 0 ) THEN 661 637 surfaces%npoints = surfaces%npoints + 1 662 point_index(k+1,j,i) = point_index_count 663 point_index_count = point_index_count + 1 664 surfaces%points(1,surfaces%npoints) = ( i  0.5_wp ) * dx 665 surfaces%points(2,surfaces%npoints) = ( j  0.5_wp ) * dy 638 point_index(k,j,i+1) = point_index_count 639 point_index_count = point_index_count + 1 640 surfaces%points(1,surfaces%npoints) = ( i + 1  0.5_wp ) * dx 641 surfaces%points(2,surfaces%npoints) = ( j  0.5_wp ) * dy 642 surfaces%points(3,surfaces%npoints) = zw(k) 643 ENDIF 644 IF ( point_index(k,j+1,i+1) < 0 ) THEN 645 surfaces%npoints = surfaces%npoints + 1 646 point_index(k,j+1,i+1) = point_index_count 647 point_index_count = point_index_count + 1 648 surfaces%points(1,surfaces%npoints) = ( i + 1  0.5_wp ) * dx 649 surfaces%points(2,surfaces%npoints) = ( j + 1  0.5_wp ) * dy 650 surfaces%points(3,surfaces%npoints) = zw(k) 651 ENDIF 652 IF ( point_index(k,j+1,i) < 0 ) THEN 653 surfaces%npoints = surfaces%npoints + 1 654 point_index(k,j+1,i) = point_index_count 655 point_index_count = point_index_count + 1 656 surfaces%points(1,surfaces%npoints) = ( i  0.5_wp ) * dx 657 surfaces%points(2,surfaces%npoints) = ( j + 1  0.5_wp ) * dy 666 658 surfaces%points(3,surfaces%npoints) = zw(k) 667 659 ENDIF 668 660 669 npg = npg + 1 670 IF ( l == 0 .OR. l == 1 ) THEN 671 surfaces%polygons(1,npg) = 4 672 surfaces%polygons(2,npg) = point_index(k,j,i) 673 surfaces%polygons(3,npg) = point_index(k,j,i+1) 674 surfaces%polygons(4,npg) = point_index(k+1,j,i+1) 675 surfaces%polygons(5,npg) = point_index(k+1,j,i) 676 ELSE 677 surfaces%polygons(1,npg) = 4 678 surfaces%polygons(2,npg) = point_index(k,j,i) 679 surfaces%polygons(3,npg) = point_index(k,j+1,i) 680 surfaces%polygons(4,npg) = point_index(k+1,j+1,i) 681 surfaces%polygons(5,npg) = point_index(k+1,j,i) 682 ENDIF 683 661 npg = npg + 1 662 surfaces%polygons(1,npg) = 4 663 surfaces%polygons(2,npg) = point_index(k,j,i) 664 surfaces%polygons(3,npg) = point_index(k,j,i+1) 665 surfaces%polygons(4,npg) = point_index(k,j+1,i+1) 666 surfaces%polygons(5,npg) = point_index(k,j+1,i) 684 667 ENDDO 685 668 686 DO m = 1, surf_lsm_v(l)%ns 687 i = surf_lsm_v(l)%i(m) + MERGE( surf_lsm_v(l)%ioff, 0, l == 3 ) 688 j = surf_lsm_v(l)%j(m) + MERGE( surf_lsm_v(l)%joff, 0, l == 1 ) 689 k = surf_lsm_v(l)%k(m) + surf_lsm_v(l)%koff 690 ! 691 ! Lower left /front vertex 692 IF ( point_index(k,j,i) < 0 ) THEN 693 surfaces%npoints = surfaces%npoints + 1 694 point_index(k,j,i) = point_index_count 695 point_index_count = point_index_count + 1 696 surfaces%points(1,surfaces%npoints) = ( i  0.5_wp ) * dx 697 surfaces%points(2,surfaces%npoints) = ( j  0.5_wp ) * dy 698 surfaces%points(3,surfaces%npoints) = zw(k1) 699 ENDIF 700 ! 701 ! Upper / lower right index for north and southfacing surfaces 702 IF ( l == 0 .OR. l == 1 ) THEN 703 IF ( point_index(k,j,i+1) < 0 ) THEN 704 surfaces%npoints = surfaces%npoints + 1 705 point_index(k,j,i+1) = point_index_count 706 point_index_count = point_index_count + 1 707 surfaces%points(1,surfaces%npoints) = ( i + 1  0.5_wp ) * dx 708 surfaces%points(2,surfaces%npoints) = ( j  0.5_wp ) * dy 669 DO l = 0, 3 670 DO m = 1, surf_def_v(l)%ns 671 ! 672 ! Determine the indices of the respective grid cell inside the topography. 673 ! Please note, jindex for northfacing surfaces ( l==0 ) is 674 ! identical to the reference jindex outside the grid box. 675 ! Equivalent for eastfacing surfaces and iindex. 676 i = surf_def_v(l)%i(m) + MERGE( surf_def_v(l)%ioff, 0, l == 3 ) 677 j = surf_def_v(l)%j(m) + MERGE( surf_def_v(l)%joff, 0, l == 1 ) 678 k = surf_def_v(l)%k(m) + surf_def_v(l)%koff 679 ! 680 ! Lower left /front vertex 681 IF ( point_index(k,j,i) < 0 ) THEN 682 surfaces%npoints = surfaces%npoints + 1 683 point_index(k,j,i) = point_index_count 684 point_index_count = point_index_count + 1 685 surfaces%points(1,surfaces%npoints) = ( i  0.5_wp ) * dx 686 surfaces%points(2,surfaces%npoints) = ( j  0.5_wp ) * dy 709 687 surfaces%points(3,surfaces%npoints) = zw(k1) 710 688 ENDIF 711 IF ( point_index(k+1,j,i+1) < 0 ) THEN 712 surfaces%npoints = surfaces%npoints + 1 713 point_index(k+1,j,i+1) = point_index_count 714 point_index_count = point_index_count + 1 715 surfaces%points(1,surfaces%npoints) = ( i + 1  0.5_wp ) * dx 716 surfaces%points(2,surfaces%npoints) = ( j  0.5_wp ) * dy 689 ! 690 ! Upper / lower right index for north and southfacing surfaces 691 IF ( l == 0 .OR. l == 1 ) THEN 692 IF ( point_index(k,j,i+1) < 0 ) THEN 693 surfaces%npoints = surfaces%npoints + 1 694 point_index(k,j,i+1) = point_index_count 695 point_index_count = point_index_count + 1 696 surfaces%points(1,surfaces%npoints) = ( i + 1  0.5_wp ) * dx 697 surfaces%points(2,surfaces%npoints) = ( j  0.5_wp ) * dy 698 surfaces%points(3,surfaces%npoints) = zw(k1) 699 ENDIF 700 IF ( point_index(k+1,j,i+1) < 0 ) THEN 701 surfaces%npoints = surfaces%npoints + 1 702 point_index(k+1,j,i+1) = point_index_count 703 point_index_count = point_index_count + 1 704 surfaces%points(1,surfaces%npoints) = ( i + 1  0.5_wp ) * dx 705 surfaces%points(2,surfaces%npoints) = ( j  0.5_wp ) * dy 706 surfaces%points(3,surfaces%npoints) = zw(k) 707 ENDIF 708 ! 709 ! Upper / lower front index for east and westfacing surfaces 710 ELSEIF ( l == 2 .OR. l == 3 ) THEN 711 IF ( point_index(k,j+1,i) < 0 ) THEN 712 surfaces%npoints = surfaces%npoints + 1 713 point_index(k,j+1,i) = point_index_count 714 point_index_count = point_index_count + 1 715 surfaces%points(1,surfaces%npoints) = ( i  0.5_wp ) * dx 716 surfaces%points(2,surfaces%npoints) = ( j + 1  0.5_wp ) * dy 717 surfaces%points(3,surfaces%npoints) = zw(k1) 718 ENDIF 719 IF ( point_index(k+1,j+1,i) < 0 ) THEN 720 surfaces%npoints = surfaces%npoints + 1 721 point_index(k+1,j+1,i) = point_index_count 722 point_index_count = point_index_count + 1 723 surfaces%points(1,surfaces%npoints) = ( i  0.5_wp ) * dx 724 surfaces%points(2,surfaces%npoints) = ( j + 1  0.5_wp ) * dy 725 surfaces%points(3,surfaces%npoints) = zw(k) 726 ENDIF 727 ENDIF 728 ! 729 ! Upper left / front vertex 730 IF ( point_index(k+1,j,i) < 0 ) THEN 731 surfaces%npoints = surfaces%npoints + 1 732 point_index(k+1,j,i) = point_index_count 733 point_index_count = point_index_count + 1 734 surfaces%points(1,surfaces%npoints) = ( i  0.5_wp ) * dx 735 surfaces%points(2,surfaces%npoints) = ( j  0.5_wp ) * dy 717 736 surfaces%points(3,surfaces%npoints) = zw(k) 718 737 ENDIF 719 ! 720 ! Upper / lower front index for east and westfacing surfaces 721 ELSEIF ( l == 2 .OR. l == 3 ) THEN 722 IF ( point_index(k,j+1,i) < 0 ) THEN 723 surfaces%npoints = surfaces%npoints + 1 724 point_index(k,j+1,i) = point_index_count 725 point_index_count = point_index_count + 1 726 surfaces%points(1,surfaces%npoints) = ( i  0.5_wp ) * dx 727 surfaces%points(2,surfaces%npoints) = ( j + 1  0.5_wp ) * dy 738 739 npg = npg + 1 740 IF ( l == 0 .OR. l == 1 ) THEN 741 surfaces%polygons(1,npg) = 4 742 surfaces%polygons(2,npg) = point_index(k,j,i) 743 surfaces%polygons(3,npg) = point_index(k,j,i+1) 744 surfaces%polygons(4,npg) = point_index(k+1,j,i+1) 745 surfaces%polygons(5,npg) = point_index(k+1,j,i) 746 ELSE 747 surfaces%polygons(1,npg) = 4 748 surfaces%polygons(2,npg) = point_index(k,j,i) 749 surfaces%polygons(3,npg) = point_index(k,j+1,i) 750 surfaces%polygons(4,npg) = point_index(k+1,j+1,i) 751 surfaces%polygons(5,npg) = point_index(k+1,j,i) 752 ENDIF 753 754 ENDDO 755 756 DO m = 1, surf_lsm_v(l)%ns 757 i = surf_lsm_v(l)%i(m) + MERGE( surf_lsm_v(l)%ioff, 0, l == 3 ) 758 j = surf_lsm_v(l)%j(m) + MERGE( surf_lsm_v(l)%joff, 0, l == 1 ) 759 k = surf_lsm_v(l)%k(m) + surf_lsm_v(l)%koff 760 ! 761 ! Lower left /front vertex 762 IF ( point_index(k,j,i) < 0 ) THEN 763 surfaces%npoints = surfaces%npoints + 1 764 point_index(k,j,i) = point_index_count 765 point_index_count = point_index_count + 1 766 surfaces%points(1,surfaces%npoints) = ( i  0.5_wp ) * dx 767 surfaces%points(2,surfaces%npoints) = ( j  0.5_wp ) * dy 728 768 surfaces%points(3,surfaces%npoints) = zw(k1) 729 769 ENDIF 730 IF ( point_index(k+1,j+1,i) < 0 ) THEN 731 surfaces%npoints = surfaces%npoints + 1 732 point_index(k+1,j+1,i) = point_index_count 733 point_index_count = point_index_count + 1 734 surfaces%points(1,surfaces%npoints) = ( i  0.5_wp ) * dx 735 surfaces%points(2,surfaces%npoints) = ( j + 1  0.5_wp ) * dy 770 ! 771 ! Upper / lower right index for north and southfacing surfaces 772 IF ( l == 0 .OR. l == 1 ) THEN 773 IF ( point_index(k,j,i+1) < 0 ) THEN 774 surfaces%npoints = surfaces%npoints + 1 775 point_index(k,j,i+1) = point_index_count 776 point_index_count = point_index_count + 1 777 surfaces%points(1,surfaces%npoints) = ( i + 1  0.5_wp ) * dx 778 surfaces%points(2,surfaces%npoints) = ( j  0.5_wp ) * dy 779 surfaces%points(3,surfaces%npoints) = zw(k1) 780 ENDIF 781 IF ( point_index(k+1,j,i+1) < 0 ) THEN 782 surfaces%npoints = surfaces%npoints + 1 783 point_index(k+1,j,i+1) = point_index_count 784 point_index_count = point_index_count + 1 785 surfaces%points(1,surfaces%npoints) = ( i + 1  0.5_wp ) * dx 786 surfaces%points(2,surfaces%npoints) = ( j  0.5_wp ) * dy 787 surfaces%points(3,surfaces%npoints) = zw(k) 788 ENDIF 789 ! 790 ! Upper / lower front index for east and westfacing surfaces 791 ELSEIF ( l == 2 .OR. l == 3 ) THEN 792 IF ( point_index(k,j+1,i) < 0 ) THEN 793 surfaces%npoints = surfaces%npoints + 1 794 point_index(k,j+1,i) = point_index_count 795 point_index_count = point_index_count + 1 796 surfaces%points(1,surfaces%npoints) = ( i  0.5_wp ) * dx 797 surfaces%points(2,surfaces%npoints) = ( j + 1  0.5_wp ) * dy 798 surfaces%points(3,surfaces%npoints) = zw(k1) 799 ENDIF 800 IF ( point_index(k+1,j+1,i) < 0 ) THEN 801 surfaces%npoints = surfaces%npoints + 1 802 point_index(k+1,j+1,i) = point_index_count 803 point_index_count = point_index_count + 1 804 surfaces%points(1,surfaces%npoints) = ( i  0.5_wp ) * dx 805 surfaces%points(2,surfaces%npoints) = ( j + 1  0.5_wp ) * dy 806 surfaces%points(3,surfaces%npoints) = zw(k) 807 ENDIF 808 ENDIF 809 ! 810 ! Upper left / front vertex 811 IF ( point_index(k+1,j,i) < 0 ) THEN 812 surfaces%npoints = surfaces%npoints + 1 813 point_index(k+1,j,i) = point_index_count 814 point_index_count = point_index_count + 1 815 surfaces%points(1,surfaces%npoints) = ( i  0.5_wp ) * dx 816 surfaces%points(2,surfaces%npoints) = ( j  0.5_wp ) * dy 736 817 surfaces%points(3,surfaces%npoints) = zw(k) 737 818 ENDIF 738 ENDIF 739 ! 740 ! Upper left / front vertex 741 IF ( point_index(k+1,j,i) < 0 ) THEN 742 surfaces%npoints = surfaces%npoints + 1 743 point_index(k+1,j,i) = point_index_count 744 point_index_count = point_index_count + 1 745 surfaces%points(1,surfaces%npoints) = ( i  0.5_wp ) * dx 746 surfaces%points(2,surfaces%npoints) = ( j  0.5_wp ) * dy 747 surfaces%points(3,surfaces%npoints) = zw(k) 748 ENDIF 749 750 npg = npg + 1 751 IF ( l == 0 .OR. l == 1 ) THEN 752 surfaces%polygons(1,npg) = 4 753 surfaces%polygons(2,npg) = point_index(k,j,i) 754 surfaces%polygons(3,npg) = point_index(k,j,i+1) 755 surfaces%polygons(4,npg) = point_index(k+1,j,i+1) 756 surfaces%polygons(5,npg) = point_index(k+1,j,i) 757 ELSE 758 surfaces%polygons(1,npg) = 4 759 surfaces%polygons(2,npg) = point_index(k,j,i) 760 surfaces%polygons(3,npg) = point_index(k,j+1,i) 761 surfaces%polygons(4,npg) = point_index(k+1,j+1,i) 762 surfaces%polygons(5,npg) = point_index(k+1,j,i) 763 ENDIF 764 ENDDO 765 DO m = 1, surf_usm_v(l)%ns 766 i = surf_usm_v(l)%i(m) + MERGE( surf_usm_v(l)%ioff, 0, l == 3 ) 767 j = surf_usm_v(l)%j(m) + MERGE( surf_usm_v(l)%joff, 0, l == 1 ) 768 k = surf_usm_v(l)%k(m) + surf_usm_v(l)%koff 769 ! 770 ! Lower left /front vertex 771 IF ( point_index(k,j,i) < 0 ) THEN 772 surfaces%npoints = surfaces%npoints + 1 773 point_index(k,j,i) = point_index_count 774 point_index_count = point_index_count + 1 775 surfaces%points(1,surfaces%npoints) = ( i  0.5_wp ) * dx 776 surfaces%points(2,surfaces%npoints) = ( j  0.5_wp ) * dy 777 surfaces%points(3,surfaces%npoints) = zw(k1) 778 ENDIF 779 ! 780 ! Upper / lower right index for north and southfacing surfaces 781 IF ( l == 0 .OR. l == 1 ) THEN 782 IF ( point_index(k,j,i+1) < 0 ) THEN 783 surfaces%npoints = surfaces%npoints + 1 784 point_index(k,j,i+1) = point_index_count 785 point_index_count = point_index_count + 1 786 surfaces%points(1,surfaces%npoints) = ( i + 1  0.5_wp ) * dx 787 surfaces%points(2,surfaces%npoints) = ( j  0.5_wp ) * dy 819 820 npg = npg + 1 821 IF ( l == 0 .OR. l == 1 ) THEN 822 surfaces%polygons(1,npg) = 4 823 surfaces%polygons(2,npg) = point_index(k,j,i) 824 surfaces%polygons(3,npg) = point_index(k,j,i+1) 825 surfaces%polygons(4,npg) = point_index(k+1,j,i+1) 826 surfaces%polygons(5,npg) = point_index(k+1,j,i) 827 ELSE 828 surfaces%polygons(1,npg) = 4 829 surfaces%polygons(2,npg) = point_index(k,j,i) 830 surfaces%polygons(3,npg) = point_index(k,j+1,i) 831 surfaces%polygons(4,npg) = point_index(k+1,j+1,i) 832 surfaces%polygons(5,npg) = point_index(k+1,j,i) 833 ENDIF 834 ENDDO 835 DO m = 1, surf_usm_v(l)%ns 836 i = surf_usm_v(l)%i(m) + MERGE( surf_usm_v(l)%ioff, 0, l == 3 ) 837 j = surf_usm_v(l)%j(m) + MERGE( surf_usm_v(l)%joff, 0, l == 1 ) 838 k = surf_usm_v(l)%k(m) + surf_usm_v(l)%koff 839 ! 840 ! Lower left /front vertex 841 IF ( point_index(k,j,i) < 0 ) THEN 842 surfaces%npoints = surfaces%npoints + 1 843 point_index(k,j,i) = point_index_count 844 point_index_count = point_index_count + 1 845 surfaces%points(1,surfaces%npoints) = ( i  0.5_wp ) * dx 846 surfaces%points(2,surfaces%npoints) = ( j  0.5_wp ) * dy 788 847 surfaces%points(3,surfaces%npoints) = zw(k1) 789 848 ENDIF 790 IF ( point_index(k+1,j,i+1) < 0 ) THEN 791 surfaces%npoints = surfaces%npoints + 1 792 point_index(k+1,j,i+1) = point_index_count 793 point_index_count = point_index_count + 1 794 surfaces%points(1,surfaces%npoints) = ( i + 1  0.5_wp ) * dx 795 surfaces%points(2,surfaces%npoints) = ( j  0.5_wp ) * dy 849 ! 850 ! Upper / lower right index for north and southfacing surfaces 851 IF ( l == 0 .OR. l == 1 ) THEN 852 IF ( point_index(k,j,i+1) < 0 ) THEN 853 surfaces%npoints = surfaces%npoints + 1 854 point_index(k,j,i+1) = point_index_count 855 point_index_count = point_index_count + 1 856 surfaces%points(1,surfaces%npoints) = ( i + 1  0.5_wp ) * dx 857 surfaces%points(2,surfaces%npoints) = ( j  0.5_wp ) * dy 858 surfaces%points(3,surfaces%npoints) = zw(k1) 859 ENDIF 860 IF ( point_index(k+1,j,i+1) < 0 ) THEN 861 surfaces%npoints = surfaces%npoints + 1 862 point_index(k+1,j,i+1) = point_index_count 863 point_index_count = point_index_count + 1 864 surfaces%points(1,surfaces%npoints) = ( i + 1  0.5_wp ) * dx 865 surfaces%points(2,surfaces%npoints) = ( j  0.5_wp ) * dy 866 surfaces%points(3,surfaces%npoints) = zw(k) 867 ENDIF 868 ! 869 ! Upper / lower front index for east and westfacing surfaces 870 ELSEIF ( l == 2 .OR. l == 3 ) THEN 871 IF ( point_index(k,j+1,i) < 0 ) THEN 872 surfaces%npoints = surfaces%npoints + 1 873 point_index(k,j+1,i) = point_index_count 874 point_index_count = point_index_count + 1 875 surfaces%points(1,surfaces%npoints) = ( i  0.5_wp ) * dx 876 surfaces%points(2,surfaces%npoints) = ( j + 1  0.5_wp ) * dy 877 surfaces%points(3,surfaces%npoints) = zw(k1) 878 ENDIF 879 IF ( point_index(k+1,j+1,i) < 0 ) THEN 880 surfaces%npoints = surfaces%npoints + 1 881 point_index(k+1,j+1,i) = point_index_count 882 point_index_count = point_index_count + 1 883 surfaces%points(1,surfaces%npoints) = ( i  0.5_wp ) * dx 884 surfaces%points(2,surfaces%npoints) = ( j + 1  0.5_wp ) * dy 885 surfaces%points(3,surfaces%npoints) = zw(k) 886 ENDIF 887 ENDIF 888 ! 889 ! Upper left / front vertex 890 IF ( point_index(k+1,j,i) < 0 ) THEN 891 surfaces%npoints = surfaces%npoints + 1 892 point_index(k+1,j,i) = point_index_count 893 point_index_count = point_index_count + 1 894 surfaces%points(1,surfaces%npoints) = ( i  0.5_wp ) * dx 895 surfaces%points(2,surfaces%npoints) = ( j  0.5_wp ) * dy 796 896 surfaces%points(3,surfaces%npoints) = zw(k) 797 897 ENDIF 798 ! 799 ! Upper / lower front index for east and westfacing surfaces 800 ELSEIF ( l == 2 .OR. l == 3 ) THEN 801 IF ( point_index(k,j+1,i) < 0 ) THEN 802 surfaces%npoints = surfaces%npoints + 1 803 point_index(k,j+1,i) = point_index_count 804 point_index_count = point_index_count + 1 805 surfaces%points(1,surfaces%npoints) = ( i  0.5_wp ) * dx 806 surfaces%points(2,surfaces%npoints) = ( j + 1  0.5_wp ) * dy 807 surfaces%points(3,surfaces%npoints) = zw(k1) 808 ENDIF 809 IF ( point_index(k+1,j+1,i) < 0 ) THEN 810 surfaces%npoints = surfaces%npoints + 1 811 point_index(k+1,j+1,i) = point_index_count 812 point_index_count = point_index_count + 1 813 surfaces%points(1,surfaces%npoints) = ( i  0.5_wp ) * dx 814 surfaces%points(2,surfaces%npoints) = ( j + 1  0.5_wp ) * dy 815 surfaces%points(3,surfaces%npoints) = zw(k) 816 ENDIF 817 ENDIF 818 ! 819 ! Upper left / front vertex 820 IF ( point_index(k+1,j,i) < 0 ) THEN 821 surfaces%npoints = surfaces%npoints + 1 822 point_index(k+1,j,i) = point_index_count 823 point_index_count = point_index_count + 1 824 surfaces%points(1,surfaces%npoints) = ( i  0.5_wp ) * dx 825 surfaces%points(2,surfaces%npoints) = ( j  0.5_wp ) * dy 826 surfaces%points(3,surfaces%npoints) = zw(k) 827 ENDIF 828 829 npg = npg + 1 830 IF ( l == 0 .OR. l == 1 ) THEN 831 surfaces%polygons(1,npg) = 4 832 surfaces%polygons(2,npg) = point_index(k,j,i) 833 surfaces%polygons(3,npg) = point_index(k,j,i+1) 834 surfaces%polygons(4,npg) = point_index(k+1,j,i+1) 835 surfaces%polygons(5,npg) = point_index(k+1,j,i) 836 ELSE 837 surfaces%polygons(1,npg) = 4 838 surfaces%polygons(2,npg) = point_index(k,j,i) 839 surfaces%polygons(3,npg) = point_index(k,j+1,i) 840 surfaces%polygons(4,npg) = point_index(k+1,j+1,i) 841 surfaces%polygons(5,npg) = point_index(k+1,j,i) 842 ENDIF 898 899 npg = npg + 1 900 IF ( l == 0 .OR. l == 1 ) THEN 901 surfaces%polygons(1,npg) = 4 902 surfaces%polygons(2,npg) = point_index(k,j,i) 903 surfaces%polygons(3,npg) = point_index(k,j,i+1) 904 surfaces%polygons(4,npg) = point_index(k+1,j,i+1) 905 surfaces%polygons(5,npg) = point_index(k+1,j,i) 906 ELSE 907 surfaces%polygons(1,npg) = 4 908 surfaces%polygons(2,npg) = point_index(k,j,i) 909 surfaces%polygons(3,npg) = point_index(k,j+1,i) 910 surfaces%polygons(4,npg) = point_index(k+1,j+1,i) 911 surfaces%polygons(5,npg) = point_index(k+1,j,i) 912 ENDIF 913 ENDDO 914 843 915 ENDDO 844 845 ENDDO 846 ! 847 ! Deallocate temporary dummy variable 848 DEALLOCATE ( point_index ) 849 ! 850 ! Sumup total number of surface elements and vertices on domain. This 851 ! will be needed for postprocessing. 852 surfaces%npoints_total = 0 916 ! 917 ! Deallocate temporary dummy variable 918 DEALLOCATE ( point_index ) 919 ! 920 ! Sumup total number of vertices on domain. This 921 ! will be needed for postprocessing. 922 surfaces%npoints_total = 0 853 923 #if defined( __parallel ) 854 CALL MPI_ALLREDUCE( surfaces%npoints, surfaces%npoints_total, 1, & 855 MPI_INTEGER, MPI_SUM, comm2d, ierr ) 856 CALL MPI_ALLREDUCE( surfaces%ns, surfaces%ns_total, 1, & 857 MPI_INTEGER, MPI_SUM, comm2d, ierr ) 924 CALL MPI_ALLREDUCE( surfaces%npoints, surfaces%npoints_total, 1, & 925 MPI_INTEGER, MPI_SUM, comm2d, ierr ) 858 926 #else 859 surfaces%npoints_total = surfaces%npoints 860 surfaces%ns_total = surfaces%ns 927 surfaces%npoints_total = surfaces%npoints 861 928 #endif 929 ENDIF 930 ! 931 ! If output to netcdf is enabled, setup the coordinate arrays that 932 ! unambiguously describe the position and orientation of each surface 933 ! element. 934 IF ( to_netcdf ) THEN 935 ! 936 ! Allocate local coordinate arrays 937 ALLOCATE( surfaces%s(1:surfaces%ns) ) 938 ALLOCATE( surfaces%xs(1:surfaces%ns) ) 939 ALLOCATE( surfaces%ys(1:surfaces%ns) ) 940 ALLOCATE( surfaces%zs(1:surfaces%ns) ) 941 ALLOCATE( surfaces%azimuth(1:surfaces%ns) ) 942 ALLOCATE( surfaces%zenith(1:surfaces%ns) ) 943 ALLOCATE( surfaces%es_utm(1:surfaces%ns) ) 944 ALLOCATE( surfaces%ns_utm(1:surfaces%ns) ) 945 ! 946 ! Gather the number of surface on each processor, in order to number 947 ! the surface elements in ascending order with respect to the total 948 ! number of surfaces in the domain. 949 #if defined( __parallel ) 950 CALL MPI_ALLGATHER( surfaces%ns, 1, MPI_INTEGER, & 951 num_surfaces_on_pe, 1, MPI_INTEGER, comm2d, ierr ) 952 #else 953 num_surfaces_on_pe = surfaces%ns 954 #endif 955 ! 956 ! First, however, determine the offset where couting of the surfaces 957 ! should start (the sum of surfaces on all PE's with lower MPI rank). 958 i = 0 959 start_count = 1 960 DO WHILE ( i < myid .AND. i <= SIZE( num_surfaces_on_pe ) ) 961 start_count = start_count + num_surfaces_on_pe(i) 962 i = i + 1 963 ENDDO 964 ! 965 ! Set coordinate arrays. For horizontal surfaces, azimuth 966 ! angles are not defined (fill value). Zenith angle is 0 (180) for 967 ! upward (downward)facing surfaces. 968 i = start_count 969 mm = 1 970 DO m = 1, surf_def_h(0)%ns 971 surfaces%s(mm) = i 972 surfaces%xs(mm) = ( surf_def_h(0)%i(m) + 0.5_wp ) * dx 973 surfaces%ys(mm) = ( surf_def_h(0)%j(m) + 0.5_wp ) * dy 974 surfaces%zs(mm) = zw(surf_def_h(0)%k(m)+surf_def_h(0)%koff) 975 surfaces%azimuth(mm) = surfaces%fillvalue 976 surfaces%zenith(mm) = 0.0 977 i = i + 1 978 mm = mm + 1 979 ENDDO 980 DO m = 1, surf_lsm_h%ns 981 surfaces%s(mm) = i 982 surfaces%xs(mm) = ( surf_lsm_h%i(m) + 0.5_wp ) * dx 983 surfaces%ys(mm) = ( surf_lsm_h%j(m) + 0.5_wp ) * dy 984 surfaces%zs(mm) = zw(surf_lsm_h%k(m)+surf_lsm_h%koff) 985 surfaces%azimuth(mm) = surfaces%fillvalue 986 surfaces%zenith(mm) = 0.0 987 i = i + 1 988 mm = mm + 1 989 ENDDO 990 DO m = 1, surf_usm_h%ns 991 surfaces%s(mm) = i 992 surfaces%xs(mm) = ( surf_usm_h%i(m) + 0.5_wp ) * dx 993 surfaces%ys(mm) = ( surf_usm_h%j(m) + 0.5_wp ) * dy 994 surfaces%zs(mm) = zw(surf_usm_h%k(m)+surf_usm_h%koff) 995 surfaces%azimuth(mm) = surfaces%fillvalue 996 surfaces%zenith(mm) = 0.0 997 i = i + 1 998 mm = mm + 1 999 ENDDO 1000 DO m = 1, surf_def_h(1)%ns 1001 surfaces%s(mm) = i 1002 surfaces%xs(mm) = ( surf_def_h(1)%i(m) + 0.5_wp ) * dx 1003 surfaces%ys(mm) = ( surf_def_h(1)%j(m) + 0.5_wp ) * dy 1004 surfaces%zs(mm) = zw(surf_def_h(1)%k(m)+surf_def_h(1)%koff) 1005 surfaces%azimuth(mm) = surfaces%fillvalue 1006 surfaces%zenith(mm) = 180.0 1007 i = i + 1 1008 mm = mm + 1 1009 ENDDO 1010 ! 1011 ! For vertical surfaces, zenith angles are not defined (fill value). 1012 ! Azimuth angle: eastward (0), westward (180), northward (270), 1013 ! southward (90). 1014 DO l = 0, 3 1015 IF ( l == 0 ) THEN 1016 az = 270.0_wp 1017 off_x = 0.0_wp 1018 off_y =  0.5_wp 1019 ELSEIF ( l == 1 ) THEN 1020 az = 90.0_wp 1021 off_x = 0.0_wp 1022 off_y = 0.5_wp 1023 ELSEIF ( l == 2 ) THEN 1024 az = 0.0_wp 1025 off_x =  0.5_wp 1026 off_y = 0.0_wp 1027 ELSEIF ( l == 3 ) THEN 1028 az = 180.0_wp 1029 off_x = 0.5_wp 1030 off_y = 0.0_wp 1031 ENDIF 1032 1033 DO m = 1, surf_def_v(l)%ns 1034 surfaces%s(mm) = i 1035 surfaces%xs(mm) = ( surf_def_v(l)%i(m) + off_x ) * dx 1036 surfaces%ys(mm) = ( surf_def_v(l)%j(m) + off_y ) * dy 1037 surfaces%zs(mm) = zu(surf_def_v(l)%k(m)) 1038 surfaces%azimuth(mm) = az 1039 surfaces%zenith(mm) = surfaces%fillvalue 1040 i = i + 1 1041 mm = mm + 1 1042 ENDDO 1043 DO m = 1, surf_lsm_v(l)%ns 1044 surfaces%s(mm) = i 1045 surfaces%xs(mm) = ( surf_lsm_v(l)%i(m) + off_x ) * dx 1046 surfaces%ys(mm) = ( surf_lsm_v(l)%j(m) + off_y ) * dy 1047 surfaces%zs(mm) = zu(surf_lsm_v(l)%k(m)) 1048 surfaces%azimuth(mm) = az 1049 surfaces%zenith(mm) = surfaces%fillvalue 1050 i = i + 1 1051 mm = mm + 1 1052 ENDDO 1053 DO m = 1, surf_usm_v(l)%ns 1054 surfaces%s(mm) = i 1055 surfaces%xs(mm) = ( surf_usm_v(l)%i(m) + off_x ) * dx 1056 surfaces%ys(mm) = ( surf_usm_v(l)%j(m) + off_y ) * dy 1057 surfaces%zs(mm) = zu(surf_usm_v(l)%k(m)) 1058 surfaces%azimuth(mm) = az 1059 surfaces%zenith(mm) = surfaces%fillvalue 1060 i = i + 1 1061 mm = mm + 1 1062 ENDDO 1063 ENDDO 1064 ! 1065 ! Finally, define UTM coordinates, which are the x/ycoordinates 1066 ! plus the origin (lowerleft coordinate of the model domain). 1067 surfaces%es_utm = surfaces%xs + init_model%origin_x 1068 surfaces%ns_utm = surfaces%ys + init_model%origin_y 1069 ! 1070 ! Initialize NetCDF data output. Please note, local start position for 1071 ! the surface elements in the NetCDF file is surfaces%s(1), while 1072 ! the number of surfaces on the subdomain is given by surfaces%ns. 1073 #if defined( __netcdf4_parallel ) 1074 1075 ! 1076 ! Calculate number of time steps to be output 1077 ntdim_surf(0) = dosurf_time_count(0) + CEILING( & 1078 ( end_time  MAX( & 1079 MERGE(skip_time_dosurf, skip_time_dosurf + spinup_time, & 1080 data_output_during_spinup ), & 1081 simulated_time_at_begin ) & 1082 ) / dt_dosurf ) 1083 1084 ntdim_surf(1) = dosurf_time_count(1) + CEILING( & 1085 ( end_time  MAX( & 1086 MERGE( skip_time_dosurf_av, skip_time_dosurf_av & 1087 + spinup_time, data_output_during_spinup ), & 1088 simulated_time_at_begin ) & 1089 ) / dt_dosurf_av ) 1090 1091 ! 1092 ! Create NetCDF4 files for parallel writing 1093 DO av = 0, 1 1094 ! 1095 ! If there is no instantaneous data (av=0) or averaged data (av=1) 1096 ! requested, do not create the corresponding NetCDF file 1097 IF ( dosurf_no(av) == 0 ) CYCLE 1098 1099 IF ( av == 0 ) THEN 1100 filename = 'SURFACE_DATA_NETCDF' // TRIM( coupling_char ) 1101 ELSE 1102 filename = 'SURFACE_DATA_AV_NETCDF' // TRIM( coupling_char ) 1103 ENDIF 1104 ! 1105 ! Open file using netCDF4/HDF5 format, parallel 1106 nc_stat = NF90_CREATE( TRIM(filename), & 1107 IOR( NF90_NOCLOBBER, & 1108 IOR( NF90_NETCDF4, NF90_MPIIO ) ), & 1109 id_set_surf(av), COMM = comm2d, INFO = MPI_INFO_NULL ) 1110 CALL netcdf_handle_error( 'surface_data_output_mod', 5550 ) 1111 1112 ! Write some global attributes 1113 IF ( av == 0 ) THEN 1114 CALL netcdf_create_global_atts( id_set_surf(av), 'surfacedata', TRIM( run_description_header ), 5551 ) 1115 time_average_text = ' ' 1116 ELSE 1117 CALL netcdf_create_global_atts( id_set_surf(av), 'surfacedata_av', TRIM( run_description_header ), 5552 ) 1118 WRITE ( time_average_text,'(F7.1,'' s avg'')' ) averaging_interval_surf 1119 nc_stat = NF90_PUT_ATT( id_set_surf(av), NF90_GLOBAL, 'time_avg', & 1120 TRIM( time_average_text ) ) 1121 CALL netcdf_handle_error( 'surface_data_output_mod', 5553 ) 1122 ENDIF 1123 1124 1125 ! 1126 ! Define time coordinate for surface data. 1127 ! For parallel output the time dimension has to be limited (ntdim_surf), 1128 ! otherwise the performance drops significantly. 1129 CALL netcdf_create_dim( id_set_surf(av), 'time', ntdim_surf(av), & 1130 id_dim_time_surf(av), 5554 ) 1131 1132 CALL netcdf_create_var( id_set_surf(av), (/ id_dim_time_surf(av) /), & 1133 'time', NF90_DOUBLE, id_var_time_surf(av), & 1134 'seconds since '//TRIM(init_model%origin_time), & 1135 'time', 5555, 5555, 5555 ) 1136 CALL netcdf_create_att( id_set_surf(av), id_var_time_surf(av), 'standard_name', 'time', 5556) 1137 CALL netcdf_create_att( id_set_surf(av), id_var_time_surf(av), 'axis', 'T', 5557) 1138 ! 1139 ! Define spatial dimensions and coordinates: 1140 ! Define index of surface element 1141 CALL netcdf_create_dim( id_set_surf(av), 's', surfaces%ns_total, & 1142 id_dim_s_surf(av), 5558 ) 1143 CALL netcdf_create_var( id_set_surf(av), (/ id_dim_s_surf(av) /), & 1144 's', NF90_DOUBLE, id_var_s_surf(av), & 1145 '1', '', 5559, 5559, 5559 ) 1146 ! 1147 ! Define x coordinate 1148 CALL netcdf_create_var( id_set_surf(av), (/ id_dim_s_surf(av) /), & 1149 'xs', NF90_DOUBLE, id_var_xs_surf(av), & 1150 'meters', '', 5561, 5561, 5561 ) 1151 ! 1152 ! Define y coordinate 1153 CALL netcdf_create_var( id_set_surf(av), (/ id_dim_s_surf(av) /), & 1154 'ys', NF90_DOUBLE, id_var_ys_surf(av), & 1155 'meters', '', 5562, 5562, 5562 ) 1156 ! 1157 ! Define z coordinate 1158 CALL netcdf_create_var( id_set_surf(av), (/ id_dim_s_surf(av) /), & 1159 'zs', NF90_DOUBLE, id_var_zs_surf(av), & 1160 'meters', '', 5560, 5560, 5560 ) 1161 1162 ! 1163 ! Define UTM coordinates 1164 CALL netcdf_create_var( id_set_surf(av), (/ id_dim_s_surf(av) /), & 1165 'Es_UTM', NF90_DOUBLE, id_var_etum_surf(av), & 1166 'meters', '', 5563, 5563, 5563 ) 1167 CALL netcdf_create_var( id_set_surf(av), (/ id_dim_s_surf(av) /), & 1168 'Ns_UTM', NF90_DOUBLE, id_var_nutm_surf(av), & 1169 'meters', '', 5564, 5564, 5564 ) 1170 ! 1171 ! Define the variables 1172 var_list = ';' 1173 i = 1 1174 1175 DO WHILE ( dosurf(av,i)(1:1) /= ' ' ) 1176 1177 CALL netcdf_create_var( id_set_surf(av),(/ id_dim_s_surf(av), & 1178 id_dim_time_surf(av) /), dosurf(av,i), & 1179 NF90_REAL4, id_var_dosurf(av,i), & 1180 dosurf_unit(av,i), dosurf(av,i), 5565, & 1181 5565, 5565, .TRUE. ) 1182 ! 1183 ! Set no fill for every variable to increase performance. 1184 nc_stat = NF90_DEF_VAR_FILL( id_set_surf(av), & 1185 id_var_dosurf(av,i), & 1186 1, 0 ) 1187 CALL netcdf_handle_error( 'surface_data_output_init', 5566 ) 1188 ! 1189 ! Set collective io operations for parallel io 1190 nc_stat = NF90_VAR_PAR_ACCESS( id_set_surf(av), & 1191 id_var_dosurf(av,i), & 1192 NF90_COLLECTIVE ) 1193 CALL netcdf_handle_error( 'surface_data_output_init', 5567 ) 1194 var_list = TRIM( var_list ) // TRIM( dosurf(av,i) ) // ';' 1195 1196 i = i + 1 1197 1198 ENDDO 1199 ! 1200 ! Write the list of variables as global attribute (this is used by 1201 ! restart runs and by combine_plot_fields) 1202 nc_stat = NF90_PUT_ATT( id_set_surf(av), NF90_GLOBAL, 'VAR_LIST', & 1203 var_list ) 1204 CALL netcdf_handle_error( 'surface_data_output_init', 5568 ) 1205 1206 ! 1207 ! Set general no fill, otherwise the performance drops significantly for 1208 ! parallel output. 1209 nc_stat = NF90_SET_FILL( id_set_surf(av), NF90_NOFILL, oldmode ) 1210 CALL netcdf_handle_error( 'surface_data_output_init', 5569 ) 1211 1212 ! 1213 ! Leave netCDF define mode 1214 nc_stat = NF90_ENDDEF( id_set_surf(av) ) 1215 CALL netcdf_handle_error( 'surface_data_output_init', 5570 ) 1216 1217 ! 1218 ! These data are only written by PE0 for parallel output to increase 1219 ! the performance. 1220 IF ( myid == 0 ) THEN 1221 ! 1222 ! Write data for surface indices 1223 ALLOCATE( netcdf_data_1d(1:surfaces%ns_total) ) 1224 1225 DO i = 1, surfaces%ns_total 1226 netcdf_data_1d(i) = i 1227 ENDDO 1228 1229 nc_stat = NF90_PUT_VAR( id_set_surf(av), id_var_s_surf(av), & 1230 netcdf_data_1d, start = (/ 1 /), & 1231 count = (/ surfaces%ns_total /) ) 1232 CALL netcdf_handle_error( 'surface_data_output_init', 5571 ) 1233 1234 DEALLOCATE( netcdf_data_1d ) 1235 1236 ENDIF 1237 1238 ! 1239 ! Write surface positions to file 1240 nc_stat = NF90_PUT_VAR( id_set_surf(av), id_var_xs_surf(av), & 1241 surfaces%xs, start = (/ surfaces%s(1) /), & 1242 count = (/ surfaces%ns /) ) 1243 CALL netcdf_handle_error( 'surface_data_output_init', 5572 ) 1244 1245 nc_stat = NF90_PUT_VAR( id_set_surf(av), id_var_ys_surf(av), & 1246 surfaces%ys, start = (/ surfaces%s(1) /), & 1247 count = (/ surfaces%ns /) ) 1248 CALL netcdf_handle_error( 'surface_data_output_init', 5573 ) 1249 1250 nc_stat = NF90_PUT_VAR( id_set_surf(av), id_var_zs_surf(av), & 1251 surfaces%zs, start = (/ surfaces%s(1) /), & 1252 count = (/ surfaces%ns /) ) 1253 CALL netcdf_handle_error( 'surface_data_output_init', 5574 ) 1254 1255 nc_stat = NF90_PUT_VAR( id_set_surf(av), id_var_etum_surf(av), & 1256 surfaces%es_utm, start = (/ surfaces%s(1) /), & 1257 count = (/ surfaces%ns /) ) 1258 CALL netcdf_handle_error( 'surface_data_output_init', 5575 ) 1259 1260 nc_stat = NF90_PUT_VAR( id_set_surf(av), id_var_nutm_surf(av), & 1261 surfaces%ns_utm, start = (/ surfaces%s(1) /), & 1262 count = (/ surfaces%ns /) ) 1263 CALL netcdf_handle_error( 'surface_data_output_init', 5576 ) 1264 1265 ENDDO 1266 #endif 1267 1268 ENDIF 862 1269 863 1270 END SUBROUTINE surface_data_output_init … … 892 1299 IF ( dosurf_no(av) == 0 ) RETURN 893 1300 ! 894 ! Open files 895 CALL check_open( 25+av ) 896 ! 897 ! Write coordinates 898 IF ( .NOT. first_output(av) ) THEN 899 DO i = 0, io_blocks1 900 IF ( i == io_group ) THEN 901 WRITE ( 25+av ) surfaces%npoints 902 WRITE ( 25+av ) surfaces%npoints_total 903 WRITE ( 25+av ) surfaces%ns 904 WRITE ( 25+av ) surfaces%ns_total 905 WRITE ( 25+av ) surfaces%points 906 WRITE ( 25+av ) surfaces%polygons 907 ENDIF 1301 ! In case of VTK output, check if binary files are open and write coordinates. 1302 IF ( to_vtk ) THEN 1303 1304 CALL check_open( 25+av ) 1305 1306 IF ( .NOT. first_output(av) ) THEN 1307 DO i = 0, io_blocks1 1308 IF ( i == io_group ) THEN 1309 WRITE ( 25+av ) surfaces%npoints 1310 WRITE ( 25+av ) surfaces%npoints_total 1311 WRITE ( 25+av ) surfaces%ns 1312 WRITE ( 25+av ) surfaces%ns_total 1313 WRITE ( 25+av ) surfaces%points 1314 WRITE ( 25+av ) surfaces%polygons 1315 ENDIF 908 1316 #if defined( __parallel ) 909 CALL MPI_BARRIER( comm2d, ierr )1317 CALL MPI_BARRIER( comm2d, ierr ) 910 1318 #endif 911 first_output(av) = .TRUE. 912 ENDDO 1319 first_output(av) = .TRUE. 1320 ENDDO 1321 ENDIF 913 1322 ENDIF 914 1323 ! 1324 ! In case of NetCDF output, check if enough time steps are available in file 1325 ! and update time axis. 1326 IF ( to_netcdf ) THEN 1327 IF ( dosurf_time_count(av) + 1 > ntdim_surf(av) ) THEN 1328 WRITE ( message_string, * ) 'Output of surface data is not given at t=', & 1329 time_since_reference_point, 's because the maximum ', & 1330 'number of output time levels is exceeded.' 1331 CALL message( 'surface_data_output', 'PA0???', 0, 1, 0, 6, 0 ) 1332 1333 RETURN 1334 1335 ENDIF 1336 ! 1337 ! Update the netCDF time axis 1338 ! In case of parallel output, this is only done by PE0 to increase the 1339 ! performance. 1340 dosurf_time_count(av) = dosurf_time_count(av) + 1 1341 IF ( myid == 0 ) THEN 1342 nc_stat = NF90_PUT_VAR( id_set_surf(av), id_var_time_surf(av), & 1343 (/ time_since_reference_point /), & 1344 start = (/ dosurf_time_count(av) /), & 1345 count = (/ 1 /) ) 1346 CALL netcdf_handle_error( 'surface_data_output', 6666 ) 1347 ENDIF 1348 ENDIF 1349 1350 ! 1351 ! Cycle through output quantities and write them to file. 915 1352 n_out = 0 916 1353 DO WHILE ( dosurf(av,n_out+1)(1:1) /= ' ' ) … … 2269 2706 !  Dimension: 1nsurfaces, 1npoints  can be ordered consecutively 2270 2707 !  Distinguish between average and nonaverage data 2271 DO i = 0, io_blocks1 2272 IF ( i == io_group ) THEN 2273 WRITE ( 25+av ) LEN_TRIM( 'time' ) 2274 WRITE ( 25+av ) 'time' 2275 WRITE ( 25+av ) time_since_reference_point 2276 WRITE ( 25+av ) LEN_TRIM( trimvar ) 2277 WRITE ( 25+av ) TRIM( trimvar ) 2278 WRITE ( 25+av ) surfaces%var_out 2279 ENDIF 2708 IF ( to_vtk ) THEN 2709 DO i = 0, io_blocks1 2710 IF ( i == io_group ) THEN 2711 WRITE ( 25+av ) LEN_TRIM( 'time' ) 2712 WRITE ( 25+av ) 'time' 2713 WRITE ( 25+av ) time_since_reference_point 2714 WRITE ( 25+av ) LEN_TRIM( trimvar ) 2715 WRITE ( 25+av ) TRIM( trimvar ) 2716 WRITE ( 25+av ) surfaces%var_out 2717 ENDIF 2280 2718 #if defined( __parallel ) 2281 CALL MPI_BARRIER( comm2d, ierr )2719 CALL MPI_BARRIER( comm2d, ierr ) 2282 2720 #endif 2283 ENDDO 2721 ENDDO 2722 ENDIF 2723 2724 IF ( to_netcdf ) THEN 2725 ! 2726 ! Write output array to file 2727 nc_stat = NF90_PUT_VAR( id_set_surf(av), id_var_dosurf(av,n_out), & 2728 surfaces%var_out, & 2729 start = (/ surfaces%s(1), dosurf_time_count(av) /), & 2730 count = (/ surfaces%ns, 1 /) ) 2731 CALL netcdf_handle_error( 'surface_data_output', 6667 ) 2732 2733 ENDIF 2284 2734 2285 2735 ENDDO … … 3553 4003 averaging_interval_surf, data_output_surf, & 3554 4004 dt_dosurf, dt_dosurf_av, & 3555 skip_time_dosurf, skip_time_dosurf_av 4005 skip_time_dosurf, skip_time_dosurf_av, & 4006 to_netcdf, to_vtk 3556 4007 3557 4008 line = ' ' … … 3599 4050 3600 4051 CHARACTER(LEN=100) :: trimvar !< dummy for single output variable 4052 CHARACTER(LEN=100) :: unit !< dummy for unit of output variable 3601 4053 3602 4054 INTEGER(iwp) :: av !< id indicating average or nonaverage data output … … 3621 4073 'PA0087', 1, 2, 0, 6, 0 ) 3622 4074 ENDIF 3623 4075 4076 #if ! defined( __netcdf4_parallel ) 4077 ! 4078 ! Surface output via NetCDF requires parallel NetCDF 4079 IF ( to_netcdf ) THEN 4080 message_string = 'to_netcdf = .True. requires parallel NetCDF' 4081 CALL message( 'surface_data_output_check_parameters', & 4082 'PA0116', 1, 2, 0, 6, 0 ) 4083 ENDIF 4084 #endif 3624 4085 ! 3625 4086 ! Count number of output variables and separate output strings for … … 3646 4107 3647 4108 ! 3648 ! Check if all output variables are known 4109 ! Check if all output variables are known and assign a unit 4110 unit = 'not set' 3649 4111 SELECT CASE ( TRIM( trimvar ) ) 3650 4112 … … 3655 4117 'PA0087', 1, 2, 0, 6, 0 ) 3656 4118 3657 CASE ( 'us', 'ts', 'qs', 'ss', 'qcs', 'ncs', 'qrs', 'nrs', 'ol' ) 3658 3659 CASE ( 'z0', 'z0h', 'z0q' ) 3660 3661 CASE ( 'theta1', 'qv1', 'thetav1' ) 4119 CASE ( 'us', 'uvw1' ) 4120 unit = 'm/s' 4121 4122 CASE ( 'ss', 'qcs', 'ncs', 'qrs', 'nrs' ) 4123 unit = '1' 4124 4125 CASE ( 'z0', 'z0h', 'z0q', 'ol' ) 4126 unit = 'm' 4127 4128 CASE ( 'ts', 'theta1', 'thetav1', 'theta_surface', 'thetav_surface' ) 4129 unit = 'K' 3662 4130 3663 4131 CASE ( 'usws', 'vsws' ) 4132 unit = 'm2/s2' 3664 4133 3665 CASE ( 'uvw1' )3666 3667 4134 CASE ( 'qcsws', 'ncsws', 'qrsws', 'nrsws', 'sasws' ) 3668 4135 3669 CASE ( 'shf', 'qsws', 'ssws' ) 3670 3671 CASE ( 'q_surface', 'theta_surface', 'thetav_surface' ) 4136 CASE ( 'shf' ) 4137 unit = 'K m/s' 4138 4139 CASE ( 'qsws' ) 4140 unit = 'kg/kg m/s' 4141 4142 CASE ( 'ssws' ) 4143 unit = 'kg/m2/s' 4144 4145 CASE ( 'qs', 'q_surface', 'qv1' ) 4146 unit = 'kg/kg' 3672 4147 3673 4148 CASE ( 'rad_net' ) 4149 unit = 'W/m2' 3674 4150 3675 4151 CASE ( 'rad_lw_in', 'rad_lw_out', 'rad_lw_dif', 'rad_lw_ref', & 3676 4152 'rad_lw_res' ) 4153 unit = 'W/m2' 3677 4154 3678 4155 CASE ( 'rad_sw_in', 'rad_sw_out', 'rad_sw_dif', 'rad_sw_ref', & 3679 4156 'rad_sw_res', 'rad_sw_dir' ) 4157 unit = 'W/m2' 3680 4158 3681 4159 CASE ( 'ghf' ) 4160 unit = 'W/m2' 3682 4161 3683 4162 CASE ( 'r_a', 'r_canopy', 'r_soil', 'r_s' ) 4163 unit = 's/m' 3684 4164 3685 4165 CASE DEFAULT … … 3689 4169 'PA0087', 1, 2, 0, 6, 0 ) 3690 4170 END SELECT 4171 4172 dosurf_unit(av,dosurf_no(av)) = unit 3691 4173 3692 4174 ENDDO … … 3717 4199 IF ( dosurf_no(av) == 0 ) RETURN 3718 4200 ! 3719 ! Open files 3720 CALL check_open( 25+av ) 3721 ! 3722 ! Write time coordinate 3723 DO i = 0, io_blocks1 3724 IF ( i == io_group ) THEN 3725 WRITE ( 25+av ) LEN_TRIM( 'END' ) 3726 WRITE ( 25+av ) 'END' 3727 ENDIF 4201 ! If output to VTK files is enabled, check if files are open and write 4202 ! an endoffile statement. 4203 IF ( to_vtk ) THEN 4204 CALL check_open( 25+av ) 4205 ! 4206 ! Write time coordinate 4207 DO i = 0, io_blocks1 4208 IF ( i == io_group ) THEN 4209 WRITE ( 25+av ) LEN_TRIM( 'END' ) 4210 WRITE ( 25+av ) 'END' 4211 ENDIF 3728 4212 #if defined( __parallel ) 3729 CALL MPI_BARRIER( comm2d, ierr )4213 CALL MPI_BARRIER( comm2d, ierr ) 3730 4214 #endif 3731 ENDDO 4215 ENDDO 4216 ENDIF 3732 4217 3733 4218 END SUBROUTINE surface_data_output_last_action
Note: See TracChangeset
for help on using the changeset viewer.