Changeset 3727 for palm/trunk/SOURCE
- Timestamp:
- Feb 8, 2019 2:52:10 PM (6 years ago)
- Location:
- palm/trunk/SOURCE
- Files:
-
- 3 edited
-
Makefile (modified) (23 diffs)
-
netcdf_interface_mod.f90 (modified) (4 diffs)
-
surface_data_output_mod.f90 (modified) (22 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/Makefile
r3700 r3727 25 25 # ----------------- 26 26 # $Id$ 27 # surface_data_output_mod depends on netcdf_interface_mod 28 # 29 # 3700 2019-01-26 17:03:42Z knoop 27 30 # Rename surface_output_mod into surface_data_output_mod 28 # 31 # 29 32 # 3637 2018-12-20 01:51:36Z knoop 30 33 # Implementation of the PALM module interface 31 # 34 # 32 35 # 3634 2018-12-18 12:31:28Z knoop 33 36 # OpenACC port for SPEC 34 # 37 # 35 38 # 3579 2018-11-29 15:32:39Z suehring 36 39 # Dependency for check_parameters on nesting_offl_mod added 37 # 40 # 38 41 # 3569 2018-11-27 17:03:40Z kanani 39 42 # dom_dwd_user, Schrempf: … … 43 46 # 3525 2018-11-14 16:06:14Z kanani 44 47 # Changes related to clean-up of biometeorology (dom_dwd_user) 45 # 48 # 46 49 # 3522 2018-11-13 12:14:36Z suehring 47 50 # Dependencies for virtual measurement module added 48 # 51 # 49 52 # 3494 2018-11-06 14:51:27Z suehring 50 53 # Surface output revised 51 # 54 # 52 55 # 3474 2018-10-30 21:07:39Z kanani 53 56 # Add virtual measurement module 54 # 57 # 55 58 # 3472 2018-10-30 20:43:50Z suehring 56 59 # Add indoor model (kanani, srissman, tlang), 57 60 # minor formatting 58 # 61 # 59 62 # 3467 2018-10-30 19:05:21Z suehring 60 63 # Implementation of a new aerosol module salsa. 61 # 62 # 64 # 65 # 63 66 # 3458 2018-10-30 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 2018-10-29 18:14:31Z kanani 71 74 # Adjustment of biometeorology dependencies 72 # 75 # 73 76 # 3436 2018-10-26 18:35:15Z gronemeier 74 77 # Add surface_mod to user_data_output_mask 75 # 78 # 76 79 # 3435 2018-10-26 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 2018-10-24 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 2018-10-19 13:09:06Z raasch 86 89 # dependencies for ocean_mod fixed 87 # 90 # 88 91 # 3355 2018-10-16 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 2018-10-15 10:38:52Z suehring 95 98 # (from branch resler) 96 99 # Add biometeorology 97 # 98 # 100 # 101 # 99 102 # 3322 2018-10-09 10:02:39Z kanani 100 103 # Formatting and cleanup 101 # 104 # 102 105 # 3298 2018-10-02 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 2018-09-24 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 2018-07-24 18:17:30Z suehring 116 119 # Bugfix, add missing dependencies for multi-agent system 117 # 120 # 118 121 # 3159 2018-07-20 11:20:01Z sward 119 122 # Added multi agent system 120 # 123 # 121 124 # 3130 2018-07-16 11:08:55Z gronemeier 122 125 # add surface_layer_fluxes_mod to turbulence_closure_mod 123 # 126 # 124 127 # 3129 2018-07-16 07:45:13Z gronemeier 125 128 # add turbulence_closure_mod to parin 126 # 129 # 127 130 # 2963 2018-04-12 14:47:44Z suehring 128 # Introduce index for vegetation/wall, pavement/green-wall and water/window 131 # Introduce index for vegetation/wall, pavement/green-wall and water/window 129 132 # surfaces, for clearer access of surface fraction, albedo, emissivity, etc. . 130 # 133 # 131 134 # 2955 2018-04-09 15:14:01Z suehring 132 135 # Add log-points to measure CPU time of NetCDF data input. 133 # 136 # 134 137 # 2938 2018-03-27 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 non-cyclic lateral boundary in nesting case 138 # 141 # 139 142 # 2936 2018-03-27 14:49:27Z suehring 140 143 # Added dependencies for parent and child synchronization 141 # 144 # 142 145 # 2921 2018-03-22 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 2018-03-21 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 2018-03-02 21:45:58Z suehring 152 155 # Changed format and enforced sorting 153 # 156 # 154 157 # 2817 2018-02-19 16:32:21Z knoop 155 158 # Preliminary gust module interface implemented 156 # 159 # 157 160 # 2802 2018-02-14 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 2018-01-30 14:12:54Z suehring 162 165 # Nesting of chemical species 163 # 166 # 164 167 # 2718 2018-01-02 08:49:38Z maronga 165 168 # Corrected "Former revisions" section 166 # 169 # 167 170 # 2697 2017-12-14 17:57:20Z kanani 168 171 # Bugfix, missing dependencies 169 # 172 # 170 173 # 2696 2017-12-14 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 2017-11-13 14:04:26Z schwenkel 185 # Added diagnostic_quantities_mod 186 # 188 # Added diagnostic_quantities_mod 189 # 187 190 # 2600 2017-11-01 14:11:20Z raasch 188 191 # comment line concerning bound checks removed 189 # 192 # 190 193 # 2599 2017-11-01 13:18:45Z hellstea 191 194 # virtual_flight_mod, synthetic_turbulence_generator_mod and … … 194 197 # 2563 2017-10-19 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 2017-10-13 18:09:32Z maronga 199 # Added date_and_time_mod 200 # 202 # Added date_and_time_mod 203 # 201 204 # 2371 2017-08-24 13:01:17Z kanani 202 205 # Corrected dependencies for vertical_nesting_mod 203 # 206 # 204 207 # 2370 2017-08-23 06:11:43Z raasch 205 208 # dependency bugfix for synthetic_turbulence_generator 206 # 209 # 207 210 # 2365 2017-08-21 14:59:59Z kanani 208 211 # Added dependencies for vertical_nesting_mod 209 # 212 # 210 213 # 2339 2017-08-07 13:55:26Z gronemeier 211 214 # corrected timestamp in header 212 # 215 # 213 216 # 2338 2017-08-07 12:15:38Z gronemeier 214 217 # Modularize 1D model 215 # 218 # 216 219 # 2320 2017-07-21 12:47:43Z suehring 217 220 # -ls_forcing nudging 218 221 # +large_scale_forcing_nudging 219 # 222 # 220 223 # 2318 2017-07-20 17:27:44Z suehring 221 224 # Add further dependencies on surface_mod 222 # 225 # 223 226 # 2317 2017-07-20 17:27:19Z suehring 224 227 # Added time_integration_spinup 225 # 228 # 226 229 # 2269 2017-06-09 11:57:32Z suehring 227 230 # Add dependency in read_3d_binary 228 # 231 # 229 232 # 2263 2017-06-08 14:59:01Z schwenkel 230 233 # Implemented splitting and merging algorithm 231 # 234 # 232 235 # 2259 2017-06-08 09:09:11Z gronemeier 233 236 # Implemented synthetic turbulence generator … … 235 238 # 2256 2017-06-07 13:58:08Z suehring 236 239 # Remove ring dependency in init_pegrid 237 # 240 # 238 241 # 2238 2017-05-31 16:49:16Z suehring 239 242 # Bugfix, further missing dependency on surface_mod 240 # 243 # 241 244 # 2237 2017-05-31 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 2017-05-30 18:08:54Z suehring 246 249 # … … 248 251 # +dependencies for surface_mod 249 252 # -wall_fluxes 250 # 253 # 251 254 # 2130 2017-01-24 16:25:39Z raasch 252 255 # dependency for timestep updated … … 254 257 # 2118 2017-01-17 16:38:49Z raasch 255 258 # -cuda_fft_interfaces_mod 256 # 259 # 257 260 # 2050 2016-11-08 15:00:55Z gronemeier 258 261 # Implement turbulent outflow method 259 # 262 # 260 263 # 2007 2016-08-24 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 2016-08-20 18:45:34Z knoop 266 269 # Bugfix: added netcdf_interface to dependency list for user_init_land_surface … … 268 271 # 1986 2016-08-10 14:07:17Z gronemeier 269 272 # POSIX-calls module added 270 # 273 # 271 274 # 1972 2016-07-26 07:52:02Z maronga 272 275 # Removed some dependencies due to further modularization of land surface model 273 # 276 # 274 277 # 1957 2016-07-07 10:43:48Z suehring 275 278 # flight module added … … 279 282 # 280 283 # 1938 2016-06-13 15:26:05Z hellstea 281 # Some dependency errors corrected 282 # 284 # Some dependency errors corrected 285 # 283 286 # 1934 2016-06-13 09:46:57Z hellstea 284 287 # poismg renamed poismg_noopt, poismg_fast_mod renamed poismg_mod 285 # 288 # 286 289 # 1914 2016-05-26 14:44:07Z witha 287 290 # Added wind_turbine_model_mod … … 294 297 # 295 298 # 1850 2016-04-08 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 2016-04-07 14:23:03Z raasch 310 313 # spectrum renamed spectra_mod, depencies for spectra changed 311 # 314 # 312 315 # 1826 2016-04-07 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 2016-04-07 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 2016-04-06 15:44:20Z maronga 321 324 # Renamed land_surface_model to land_surface_model_mod. 322 # 325 # 323 326 # 1808 2016-04-05 19:44:00Z raasch 324 327 # -local_flush, -local_getenv … … 346 349 # 347 350 # 1762 2016-02-25 12:31:13Z hellstea 348 # +pmc_interface, +pmc routines 351 # +pmc_interface, +pmc routines 349 352 # 350 353 # 1747 2016-02-08 12:25:53Z raasch … … 352 355 # 353 356 # 1691 2015-10-26 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 2015-04-30 07:05:52Z maronga 358 361 # Added user_init_radiation.f90 359 # 362 # 360 363 # 1575 2015-03-27 09:56:27Z raasch 361 364 # +poismg_fast … … 363 366 # 1551 2015-03-03 14:18:16Z maronga 364 367 # Bugfix: further adjustments for the land surface model and radiation model 365 # 368 # 366 369 # 1517 2015-01-07 19:12:25Z hoffmann 367 370 # advec_s_bc added to prognostic_equations … … 369 372 # 1500 2014-12-03 17:42:41Z maronga 370 373 # Bugfix: missing adjustments for land surface model and radiation model 371 # 374 # 372 375 # 1496 2014-12-02 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 2014-10-21 10:53:05Z kanani 377 380 # plant_canopy_model-dependency 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 2014-08-02 20:10:32Z letzel 381 384 # bugfix: cpulog added to lpm_advec 382 # 385 # 383 386 # 1404 2014-05-14 09:01:39Z keck 384 387 # bugfix: dependencies added for progress_bar … … 386 389 # 1402 2014-05-09 14:25:13Z raasch 387 390 # progress_bar added 388 # 391 # 389 392 # 1400 2014-05-09 14:03:54Z knoop 390 393 # Added new module random_generator_parallel 391 # 394 # 392 395 # 1380 2014-04-28 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 2014-04-25 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 2014-04-17 12:28:49Z keck 404 407 # bugfix: cpulog added to lpm_pack_arrays … … 406 409 # 1361 2014-04-16 15:17:48Z hoffmann 407 410 # cpulog added to microphysics 408 # 411 # 409 412 # 1359 2014-04-11 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 2019-01-26 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 2019-01-23 09:57:04Z suehring 27 30 ! Add output of surface-parallel 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 !< E-UTM coordinate for NetCDF output 124 REAL(wp), DIMENSION(:), ALLOCATABLE :: ns_utm !< E-UTM coordinate for NetCDF output 125 REAL(wp), DIMENSION(:), ALLOCATABLE :: xs !< x-coordinate for NetCDF output 126 REAL(wp), DIMENSION(:), ALLOCATABLE :: ys !< y-coordinate for NetCDF output 127 REAL(wp), DIMENSION(:), ALLOCATABLE :: zs !< z-coordinate 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 non-averaged output 136 CHARACTER(LEN=100), DIMENSION(0:1,300) :: dosurf = ' ' !< internal variable which describes the output variables 137 !< and separates averaged from non-averaged 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 set-fill-mode 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 surface-data 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 non-averaged (=0) data 180 231 INTEGER(iwp) :: i !< grid index in x-direction, also running variable for counting non-average data output 181 232 INTEGER(iwp) :: j !< grid index in y-direction, also running variable for counting average data output … … 183 234 INTEGER(iwp) :: l !< running index for surface-element 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:numprocs-1) :: num_points_on_pe !< array which contains the number of points on all mpi ranks 242 INTEGER(iwp), DIMENSION(0:numprocs-1) :: 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 x-direction between the stored grid index and the actual wall 247 REAL(wp) :: off_y !< grid offset in y-direction 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 !eastward-facing 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,nys-1:nyn+1,nxl-1: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,nys-1:nyn+1,nxl-1: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, j-index for north-facing surfaces ( l==0 ) is 298 !-- identical to the reference j-index outside the grid box. 299 !-- Equivalent for east-facing surfaces and i-index. 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 south-facing 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, j-index for north-facing surfaces 368 !-- ( l==0 ) is identical to the reference j-index outside the grid 369 !-- box. Equivalent for east-facing surfaces and i-index. 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 south-facing 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 west-facing 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 south-facing 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 west-facing 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 south-facing 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 west-facing 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 west-facing 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 south-facing 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 west-facing 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 south-facing 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 west-facing 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, j-index for north-facing surfaces ( l==0 ) is 604 !-- identical to the reference j-index outside the grid box. 605 !-- Equivalent for east-facing surfaces and i-index. 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(k-1) 618 ENDIF 619 ! 620 !-- Upper / lower right index for north- and south-facing 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(k-1) 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 west-facing 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(k-1) 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(k-1) 699 ENDIF 700 ! 701 !-- Upper / lower right index for north- and south-facing 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, j-index for north-facing surfaces ( l==0 ) is 674 !-- identical to the reference j-index outside the grid box. 675 !-- Equivalent for east-facing surfaces and i-index. 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(k-1) 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 south-facing 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(k-1) 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 west-facing 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(k-1) 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 west-facing 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(k-1) 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 south-facing 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(k-1) 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 west-facing 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(k-1) 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(k-1) 778 ENDIF 779 ! 780 !-- Upper / lower right index for north- and south-facing 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(k-1) 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 south-facing 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(k-1) 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 west-facing 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(k-1) 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 west-facing 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(k-1) 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 !-- Sum-up total number of surface elements and vertices on domain. This 851 !-- will be needed for post-processing. 852 surfaces%npoints_total = 0 916 ! 917 !-- Deallocate temporary dummy variable 918 DEALLOCATE ( point_index ) 919 ! 920 !-- Sum-up total number of vertices on domain. This 921 !-- will be needed for post-processing. 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, set-up 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/y-coordinates 1066 !-- plus the origin (lower-left 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), 'surface-data', TRIM( run_description_header ), 5551 ) 1115 time_average_text = ' ' 1116 ELSE 1117 CALL netcdf_create_global_atts( id_set_surf(av), 'surface-data_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_blocks-1 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_blocks-1 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: 1-nsurfaces, 1-npoints - can be ordered consecutively 2270 2707 !-- - Distinguish between average and non-average data 2271 DO i = 0, io_blocks-1 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_blocks-1 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 non-average 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_blocks-1 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 end-of-file statement. 4203 IF ( to_vtk ) THEN 4204 CALL check_open( 25+av ) 4205 ! 4206 !-- Write time coordinate 4207 DO i = 0, io_blocks-1 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.
