 Timestamp:
 Feb 12, 2019 12:57:51 PM (5 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

palm/trunk/SOURCE/surface_data_output_mod.f90
r3735 r3736 25 25 !  26 26 ! $Id$ 27 ! Add azimuth and zenith to output file; set longname attributes; 28 ! cleanup coding layout 29 ! 30 ! 3735 20190212 09:52:40Z dom_dwd_user 27 31 !  Split initialization into initialization of arrays and further initialization 28 32 ! in order to enable reading of restart data. 29 33 !  Consider restarts in surface data averaging. 30 34 !  Correct error message numbers 31 ! 35 ! 32 36 ! 3731 20190211 13:06:27Z suehring 33 37 ! Bugfix: add cpp options 34 ! 38 ! 35 39 ! 3727 20190208 14:52:10Z gronemeier 36 40 ! Enable NetCDF output for surface data (suehring, gronemeier) 37 ! 41 ! 38 42 ! 3691 20190123 09:57:04Z suehring 39 ! Add output of surfaceparallel flow speed 40 ! 43 ! Add output of surfaceparallel flow speed 44 ! 41 45 ! 3648 20190102 16:35:46Z suehring 42 46 ! Rename module and subroutines 43 ! 47 ! 44 48 ! 3646 20181228 17:58:49Z kanani 45 49 ! Bugfix: use time_since_reference_point instead of simulated_time (relevant 46 50 ! when using wall/soil spinup) 47 ! 51 ! 48 52 ! 3614 20181210 07:05:46Z raasch 49 53 ! unused variables removed 50 ! 54 ! 51 55 ! 3572 20181128 11:40:28Z suehring 52 56 ! Added short and longwave radiation flux arrays (e.g. diffuse, direct, 53 57 ! reflected, resedual) for all surfaces (M. Salim) 54 ! 58 ! 55 59 ! 3494 20181106 14:51:27Z suehring 56 ! Bugfix in gathering surface data from different types and orientation. 60 ! Bugfix in gathering surface data from different types and orientation. 57 61 ! Output of total number of surfaces and vertices added. 58 62 ! String length is output, for more convinient postprocessing. 59 63 ! Last actions added. 60 ! 64 ! 61 65 ! 3483 20181102 14:19:26Z raasch 62 66 ! bugfix: missed directives for MPI added 63 ! 67 ! 64 68 ! 3421 20181024 18:39:32Z gronemeier 65 69 ! Add output variables … … 77 81 !  78 82 !> Generate output for surface data. 79 !> 80 !> @todo Create namelist file for postprocessing tool. 83 !> 84 !> @todo Create namelist file for postprocessing tool. 81 85 !! 82 86 … … 87 91 USE arrays_3d, & 88 92 ONLY: zu, zw 89 93 90 94 USE control_parameters, & 91 95 ONLY: coupling_char, data_output_during_spinup, end_time, & 92 96 message_string, run_description_header, simulated_time_at_begin, & 93 97 spinup_time, surface_output 94 98 95 99 USE grid_variables, & 96 100 ONLY: dx,dy 97 101 98 102 USE indices, & 99 103 ONLY: nxl, nxr, nys, nyn, nzb, nzt 100 104 101 105 #if defined( __netcdf ) 102 106 USE NETCDF … … 109 113 ONLY: netcdf_create_att, netcdf_create_dim, netcdf_create_file, & 110 114 netcdf_create_global_atts, netcdf_create_var, netcdf_handle_error 111 115 112 116 USE pegrid 113 117 114 118 USE surface_mod, & 115 119 ONLY: surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, & … … 119 123 120 124 TYPE surf_out !< data structure which contains all surfaces elements of all types on subdomain 121 125 122 126 INTEGER(iwp) :: ns !< number of surface elements on subdomain 123 127 INTEGER(iwp) :: ns_total !< total number of surface elements 124 128 INTEGER(iwp) :: npoints !< number of points / vertices which define a surface element (on subdomain) 125 129 INTEGER(iwp) :: npoints_total !< total number of points / vertices which define a surface element 126 130 127 131 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: s !< coordinate for NetCDF output, number of the surface element 128 129 REAL(wp) :: fillvalue = 9999. 9_wp !< fillvalue for surface elements which are not defined130 131 REAL(wp), DIMENSION(:), ALLOCATABLE :: azimuth !< azimuth orientation coordinate for NetCDF output 132 REAL(wp), DIMENSION(:), ALLOCATABLE :: es_utm !< EUTM coordinate for NetCDF output 133 REAL(wp), DIMENSION(:), ALLOCATABLE :: ns_utm !< EUTM coordinate for NetCDF output 134 REAL(wp), DIMENSION(:), ALLOCATABLE :: xs !< xcoordinate for NetCDF output 135 REAL(wp), DIMENSION(:), ALLOCATABLE :: ys !< ycoordinate for NetCDF output 136 REAL(wp), DIMENSION(:), ALLOCATABLE :: zs !< zcoordinate for NetCDF output 137 REAL(wp), DIMENSION(:), ALLOCATABLE :: zenith !< zenith orientation coordinate for NetCDF output 132 133 REAL(wp) :: fillvalue = 9999.0_wp !< fillvalue for surface elements which are not defined 134 135 REAL(wp), DIMENSION(:), ALLOCATABLE :: azimuth !< azimuth orientation coordinate for NetCDF output 136 REAL(wp), DIMENSION(:), ALLOCATABLE :: es_utm !< EUTM coordinate for NetCDF output 137 REAL(wp), DIMENSION(:), ALLOCATABLE :: ns_utm !< EUTM coordinate for NetCDF output 138 REAL(wp), DIMENSION(:), ALLOCATABLE :: xs !< xcoordinate for NetCDF output 139 REAL(wp), DIMENSION(:), ALLOCATABLE :: ys !< ycoordinate for NetCDF output 140 REAL(wp), DIMENSION(:), ALLOCATABLE :: zs !< zcoordinate for NetCDF output 141 REAL(wp), DIMENSION(:), ALLOCATABLE :: zenith !< zenith orientation coordinate for NetCDF output 138 142 REAL(wp), DIMENSION(:), ALLOCATABLE :: var_out !< output variables 139 143 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: var_av !< variables used for averaging … … 141 145 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: polygons !< polygon data of a surface element 142 146 END TYPE surf_out 143 147 144 148 CHARACTER(LEN=100), DIMENSION(300) :: data_output_surf = ' ' !< namelist variable which describes the output variables 145 CHARACTER(LEN=100), DIMENSION(0:1,300) :: dosurf = ' ' !< internal variable which describes the output variables 149 CHARACTER(LEN=100), DIMENSION(0:1,300) :: dosurf = ' ' !< internal variable which describes the output variables 146 150 !< and separates averaged from nonaveraged output 147 151 CHARACTER(LEN=100), DIMENSION(0:1,300) :: dosurf_unit = ' ' !< internal variable which holds the unit of the given output variable 148 152 149 153 INTEGER(iwp) :: average_count_surf = 0 !< number of ensemble members used for averaging 150 154 INTEGER(iwp) :: dosurf_no(0:1) = 0 !< number of surface output quantities 151 155 152 156 INTEGER(iwp) :: nc_stat !< error code for netcdf routines 153 INTEGER(iwp) , SAVE :: oldmode!< save old setfillmode of netcdf file (not needed, but required for routine call)157 INTEGER(iwp) :: oldmode !< save old setfillmode of netcdf file (not needed, but required for routine call) 154 158 155 159 INTEGER(iwp), DIMENSION(0:1) :: id_dim_s_surf !< netcdf ID for dimension s 156 160 INTEGER(iwp), DIMENSION(0:1) :: id_dim_time_surf !< netcdf ID for dimension time 157 161 INTEGER(iwp), DIMENSION(0:1) :: id_set_surf !< netcdf ID for file 162 INTEGER(iwp), DIMENSION(0:1) :: id_var_azimuth_surf !< netcdf ID for variable azimuth 158 163 INTEGER(iwp), DIMENSION(0:1) :: id_var_etum_surf !< netcdf ID for variable Es_UTM 159 164 INTEGER(iwp), DIMENSION(0:1) :: id_var_nutm_surf !< netcdf ID for variable Ns_UTM … … 162 167 INTEGER(iwp), DIMENSION(0:1) :: id_var_xs_surf !< netcdf ID for variable xs 163 168 INTEGER(iwp), DIMENSION(0:1) :: id_var_ys_surf !< netcdf ID for variable ys 169 INTEGER(iwp), DIMENSION(0:1) :: id_var_zenith_surf !< netcdf ID for variable zenith 164 170 INTEGER(iwp), DIMENSION(0:1) :: id_var_zs_surf !< netcdf ID for variable zs 165 171 INTEGER(iwp), DIMENSION(0:1) :: dosurf_time_count = 0 !< count of output time steps 166 172 INTEGER(iwp), DIMENSION(0:1) :: ntdim_surf !< number of output time steps 167 173 168 174 INTEGER(iwp), DIMENSION(0:1,300) :: id_var_dosurf !< netcdf ID for output variables 169 175 170 176 LOGICAL :: first_output(0:1) = .FALSE. !< true if first output was already called 171 LOGICAL :: to_netcdf = .TRUE. !< flag indicating parallel NetCDF output 172 LOGICAL :: to_vtk = .TRUE. !< flag indicating binary surfacedata output that can be further 177 LOGICAL :: to_netcdf = .TRUE. !< flag indicating parallel NetCDF output 178 LOGICAL :: to_vtk = .TRUE. !< flag indicating binary surfacedata output that can be further 173 179 !< processed to VTK format 174 180 … … 182 188 183 189 TYPE(surf_out) :: surfaces !< variable which contains all required output information 184 190 185 191 SAVE 186 192 … … 190 196 MODULE PROCEDURE surface_data_output 191 197 END INTERFACE surface_data_output 192 198 193 199 INTERFACE surface_data_output_averaging 194 200 MODULE PROCEDURE surface_data_output_averaging 195 201 END INTERFACE surface_data_output_averaging 196 202 197 203 INTERFACE surface_data_output_check_parameters 198 204 MODULE PROCEDURE surface_data_output_check_parameters 199 205 END INTERFACE surface_data_output_check_parameters 200 206 201 207 INTERFACE surface_data_output_init 202 208 MODULE PROCEDURE surface_data_output_init 203 209 END INTERFACE surface_data_output_init 204 210 205 211 INTERFACE surface_data_output_init_arrays 206 212 MODULE PROCEDURE surface_data_output_init_arrays 207 213 END INTERFACE surface_data_output_init_arrays 208 214 209 215 INTERFACE surface_data_output_last_action 210 216 MODULE PROCEDURE surface_data_output_last_action 211 217 END INTERFACE surface_data_output_last_action 212 218 213 219 INTERFACE surface_data_output_parin 214 220 MODULE PROCEDURE surface_data_output_parin 215 221 END INTERFACE surface_data_output_parin 216 222 217 223 INTERFACE surface_data_output_rrd_global 218 224 MODULE PROCEDURE surface_data_output_rrd_global 219 225 END INTERFACE surface_data_output_rrd_global 220 226 221 227 INTERFACE surface_data_output_rrd_local 222 228 MODULE PROCEDURE surface_data_output_rrd_local 223 229 END INTERFACE surface_data_output_rrd_local 224 230 225 231 INTERFACE surface_data_output_wrd_global 226 232 MODULE PROCEDURE surface_data_output_wrd_global 227 233 END INTERFACE surface_data_output_wrd_global 228 234 229 235 INTERFACE surface_data_output_wrd_local 230 236 MODULE PROCEDURE surface_data_output_wrd_local … … 250 256 !  251 257 !> This routine counts the number of surfaces on each core and allocates 252 !> arrays. 258 !> arrays. 253 259 !! 254 260 SUBROUTINE surface_data_output_init_arrays 255 261 256 262 IMPLICIT NONE 257 263 … … 259 265 ! Determine the number of surface elements on subdomain 260 266 surfaces%ns = surf_def_h(0)%ns + surf_lsm_h%ns + surf_usm_h%ns & !horizontal upwardfacing 261 + surf_def_h(1)%ns & !horizontal downardfacing 267 + surf_def_h(1)%ns & !horizontal downardfacing 262 268 + surf_def_v(0)%ns + surf_lsm_v(0)%ns + surf_usm_v(0)%ns & !northwardfacing 263 + surf_def_v(1)%ns + surf_lsm_v(1)%ns + surf_usm_v(1)%ns & !southwardfacing 264 + surf_def_v(2)%ns + surf_lsm_v(2)%ns + surf_usm_v(2)%ns & !westwardfacing 265 + surf_def_v(3)%ns + surf_lsm_v(3)%ns + surf_usm_v(3)%ns !eastwardfacing 269 + surf_def_v(1)%ns + surf_lsm_v(1)%ns + surf_usm_v(1)%ns & !southwardfacing 270 + surf_def_v(2)%ns + surf_lsm_v(2)%ns + surf_usm_v(2)%ns & !westwardfacing 271 + surf_def_v(3)%ns + surf_lsm_v(3)%ns + surf_usm_v(3)%ns !eastwardfacing 266 272 ! 267 273 ! Determine the total number of surfaces in the model domain … … 277 283 surfaces%var_out = surfaces%fillvalue 278 284 ! 279 ! If there is an output of time average output variables, allocate the 285 ! If there is an output of time average output variables, allocate the 280 286 ! required array. 281 287 IF ( dosurf_no(1) > 0 ) THEN … … 283 289 surfaces%var_av = 0.0_wp 284 290 ENDIF 285 291 286 292 END SUBROUTINE surface_data_output_init_arrays 287 288 293 294 289 295 !! 290 296 ! Description: … … 294 300 !! 295 301 SUBROUTINE surface_data_output_init 296 302 297 303 IMPLICIT NONE 298 304 … … 311 317 INTEGER(iwp) :: point_index_count !< local counter variable for point index 312 318 INTEGER(iwp) :: start_count !< local start counter for the surface index 313 319 314 320 INTEGER(iwp), DIMENSION(0:numprocs1) :: num_points_on_pe !< array which contains the number of points on all mpi ranks 315 321 INTEGER(iwp), DIMENSION(0:numprocs1) :: num_surfaces_on_pe !< array which contains the number of surfaces on all mpi ranks 316 322 INTEGER(iwp), ALLOCATABLE, DIMENSION(:,:,:) :: point_index !< dummy array used to check where the reference points for surface polygons are located 317 323 318 324 REAL(wp) :: az !< azimuth angle, indicated the vertical orientation of a surface element 319 325 REAL(wp) :: off_x !< grid offset in xdirection between the stored grid index and the actual wall … … 323 329 324 330 ! 325 ! If output to VTK format is enabled, initialize point and polygon data. 326 ! In a first step, count the number of points which are defining 331 ! If output to VTK format is enabled, initialize point and polygon data. 332 ! In a first step, count the number of points which are defining 327 333 ! the surfaces and the polygons. 328 334 IF ( to_vtk ) THEN … … 340 346 k = surf_def_h(0)%k(m) + surf_def_h(0)%koff 341 347 ! 342 ! Check if the vertices that define the surface element are already 348 ! Check if the vertices that define the surface element are already 343 349 ! defined, if not, increment the counter. 344 350 IF ( point_index(k,j,i) < 0 ) THEN 345 surfaces%npoints = surfaces%npoints + 1 351 surfaces%npoints = surfaces%npoints + 1 346 352 point_index(k,j,i) = surfaces%npoints  1 347 353 ENDIF 348 354 IF ( point_index(k,j,i+1) < 0 ) THEN 349 surfaces%npoints = surfaces%npoints + 1 355 surfaces%npoints = surfaces%npoints + 1 350 356 point_index(k,j,i+1) = surfaces%npoints  1 351 357 ENDIF 352 358 IF ( point_index(k,j+1,i+1) < 0 ) THEN 353 surfaces%npoints = surfaces%npoints + 1 359 surfaces%npoints = surfaces%npoints + 1 354 360 point_index(k,j+1,i+1) = surfaces%npoints  1 355 361 ENDIF 356 362 IF ( point_index(k,j+1,i) < 0 ) THEN 357 surfaces%npoints = surfaces%npoints + 1 363 surfaces%npoints = surfaces%npoints + 1 358 364 point_index(k,j+1,i) = surfaces%npoints  1 359 ENDIF 365 ENDIF 360 366 ENDDO 361 367 ENDDO … … 366 372 367 373 IF ( point_index(k,j,i) < 0 ) THEN 368 surfaces%npoints = surfaces%npoints + 1 374 surfaces%npoints = surfaces%npoints + 1 369 375 point_index(k,j,i) = surfaces%npoints  1 370 376 ENDIF 371 377 IF ( point_index(k,j,i+1) < 0 ) THEN 372 surfaces%npoints = surfaces%npoints + 1 378 surfaces%npoints = surfaces%npoints + 1 373 379 point_index(k,j,i+1) = surfaces%npoints  1 374 380 ENDIF 375 381 IF ( point_index(k,j+1,i+1) < 0 ) THEN 376 surfaces%npoints = surfaces%npoints + 1 382 surfaces%npoints = surfaces%npoints + 1 377 383 point_index(k,j+1,i+1) = surfaces%npoints  1 378 384 ENDIF 379 385 IF ( point_index(k,j+1,i) < 0 ) THEN 380 surfaces%npoints = surfaces%npoints + 1 386 surfaces%npoints = surfaces%npoints + 1 381 387 point_index(k,j+1,i) = surfaces%npoints  1 382 ENDIF 388 ENDIF 383 389 ENDDO 384 390 DO m = 1, surf_usm_h%ns … … 388 394 389 395 IF ( point_index(k,j,i) < 0 ) THEN 390 surfaces%npoints = surfaces%npoints + 1 396 surfaces%npoints = surfaces%npoints + 1 391 397 point_index(k,j,i) = surfaces%npoints  1 392 398 ENDIF 393 399 IF ( point_index(k,j,i+1) < 0 ) THEN 394 surfaces%npoints = surfaces%npoints + 1 400 surfaces%npoints = surfaces%npoints + 1 395 401 point_index(k,j,i+1) = surfaces%npoints  1 396 402 ENDIF 397 403 IF ( point_index(k,j+1,i+1) < 0 ) THEN 398 surfaces%npoints = surfaces%npoints + 1 404 surfaces%npoints = surfaces%npoints + 1 399 405 point_index(k,j+1,i+1) = surfaces%npoints  1 400 406 ENDIF 401 407 IF ( point_index(k,j+1,i) < 0 ) THEN 402 surfaces%npoints = surfaces%npoints + 1 408 surfaces%npoints = surfaces%npoints + 1 403 409 point_index(k,j+1,i) = surfaces%npoints  1 404 ENDIF 410 ENDIF 405 411 ENDDO 406 412 ! … … 410 416 ! 411 417 ! Determine the indices of the respective grid cell inside the 412 ! topography. Please note, jindex for northfacing surfaces 413 ! ( l==0 ) is identical to the reference jindex outside the grid 418 ! topography. Please note, jindex for northfacing surfaces 419 ! ( l==0 ) is identical to the reference jindex outside the grid 414 420 ! box. Equivalent for eastfacing surfaces and iindex. 415 421 i = surf_def_v(l)%i(m) + MERGE( surf_def_v(l)%ioff, 0, l == 3 ) … … 419 425 ! Lower left /front vertex 420 426 IF ( point_index(k,j,i) < 0 ) THEN 421 surfaces%npoints = surfaces%npoints + 1 427 surfaces%npoints = surfaces%npoints + 1 422 428 point_index(k,j,i) = surfaces%npoints  1 423 ENDIF 429 ENDIF 424 430 ! 425 431 ! Upper / lower right index for north and southfacing surfaces 426 432 IF ( l == 0 .OR. l == 1 ) THEN 427 433 IF ( point_index(k,j,i+1) < 0 ) THEN 428 surfaces%npoints = surfaces%npoints + 1 434 surfaces%npoints = surfaces%npoints + 1 429 435 point_index(k,j,i+1) = surfaces%npoints  1 430 ENDIF 436 ENDIF 431 437 IF ( point_index(k+1,j,i+1) < 0 ) THEN 432 surfaces%npoints = surfaces%npoints + 1 438 surfaces%npoints = surfaces%npoints + 1 433 439 point_index(k+1,j,i+1) = surfaces%npoints  1 434 440 ENDIF … … 437 443 ELSEIF ( l == 2 .OR. l == 3 ) THEN 438 444 IF ( point_index(k,j+1,i) < 0 ) THEN 439 surfaces%npoints = surfaces%npoints + 1 445 surfaces%npoints = surfaces%npoints + 1 440 446 point_index(k,j+1,i) = surfaces%npoints  1 441 ENDIF 447 ENDIF 442 448 IF ( point_index(k+1,j+1,i) < 0 ) THEN 443 surfaces%npoints = surfaces%npoints + 1 449 surfaces%npoints = surfaces%npoints + 1 444 450 point_index(k+1,j+1,i) = surfaces%npoints  1 445 451 ENDIF … … 448 454 ! Upper left / front vertex 449 455 IF ( point_index(k+1,j,i) < 0 ) THEN 450 surfaces%npoints = surfaces%npoints + 1 456 surfaces%npoints = surfaces%npoints + 1 451 457 point_index(k+1,j,i) = surfaces%npoints  1 452 458 ENDIF … … 459 465 ! Lower left /front vertex 460 466 IF ( point_index(k,j,i) < 0 ) THEN 461 surfaces%npoints = surfaces%npoints + 1 467 surfaces%npoints = surfaces%npoints + 1 462 468 point_index(k,j,i) = surfaces%npoints  1 463 ENDIF 469 ENDIF 464 470 ! 465 471 ! Upper / lower right index for north and southfacing surfaces 466 472 IF ( l == 0 .OR. l == 1 ) THEN 467 473 IF ( point_index(k,j,i+1) < 0 ) THEN 468 surfaces%npoints = surfaces%npoints + 1 474 surfaces%npoints = surfaces%npoints + 1 469 475 point_index(k,j,i+1) = surfaces%npoints  1 470 ENDIF 476 ENDIF 471 477 IF ( point_index(k+1,j,i+1) < 0 ) THEN 472 surfaces%npoints = surfaces%npoints + 1 478 surfaces%npoints = surfaces%npoints + 1 473 479 point_index(k+1,j,i+1) = surfaces%npoints  1 474 480 ENDIF … … 477 483 ELSEIF ( l == 2 .OR. l == 3 ) THEN 478 484 IF ( point_index(k,j+1,i) < 0 ) THEN 479 surfaces%npoints = surfaces%npoints + 1 485 surfaces%npoints = surfaces%npoints + 1 480 486 point_index(k,j+1,i) = surfaces%npoints  1 481 ENDIF 487 ENDIF 482 488 IF ( point_index(k+1,j+1,i) < 0 ) THEN 483 surfaces%npoints = surfaces%npoints + 1 489 surfaces%npoints = surfaces%npoints + 1 484 490 point_index(k+1,j+1,i) = surfaces%npoints  1 485 491 ENDIF … … 488 494 ! Upper left / front vertex 489 495 IF ( point_index(k+1,j,i) < 0 ) THEN 490 surfaces%npoints = surfaces%npoints + 1 496 surfaces%npoints = surfaces%npoints + 1 491 497 point_index(k+1,j,i) = surfaces%npoints  1 492 498 ENDIF … … 500 506 ! Lower left /front vertex 501 507 IF ( point_index(k,j,i) < 0 ) THEN 502 surfaces%npoints = surfaces%npoints + 1 508 surfaces%npoints = surfaces%npoints + 1 503 509 point_index(k,j,i) = surfaces%npoints  1 504 ENDIF 510 ENDIF 505 511 ! 506 512 ! Upper / lower right index for north and southfacing surfaces 507 513 IF ( l == 0 .OR. l == 1 ) THEN 508 514 IF ( point_index(k,j,i+1) < 0 ) THEN 509 surfaces%npoints = surfaces%npoints + 1 515 surfaces%npoints = surfaces%npoints + 1 510 516 point_index(k,j,i+1) = surfaces%npoints  1 511 ENDIF 517 ENDIF 512 518 IF ( point_index(k+1,j,i+1) < 0 ) THEN 513 surfaces%npoints = surfaces%npoints + 1 519 surfaces%npoints = surfaces%npoints + 1 514 520 point_index(k+1,j,i+1) = surfaces%npoints  1 515 521 ENDIF … … 518 524 ELSEIF ( l == 2 .OR. l == 3 ) THEN 519 525 IF ( point_index(k,j+1,i) < 0 ) THEN 520 surfaces%npoints = surfaces%npoints + 1 526 surfaces%npoints = surfaces%npoints + 1 521 527 point_index(k,j+1,i) = surfaces%npoints  1 522 ENDIF 528 ENDIF 523 529 IF ( point_index(k+1,j+1,i) < 0 ) THEN 524 surfaces%npoints = surfaces%npoints + 1 530 surfaces%npoints = surfaces%npoints + 1 525 531 point_index(k+1,j+1,i) = surfaces%npoints  1 526 532 ENDIF … … 529 535 ! Upper left / front vertex 530 536 IF ( point_index(k+1,j,i) < 0 ) THEN 531 surfaces%npoints = surfaces%npoints + 1 537 surfaces%npoints = surfaces%npoints + 1 532 538 point_index(k+1,j,i) = surfaces%npoints  1 533 539 ENDIF … … 536 542 ENDDO 537 543 ! 538 ! Allocate the number of points and polygons. Note, the number of polygons539 ! is identical to the number of surfaces elements, whereas the number540 ! of points (vertices), which define the polygons, can be larger.544 ! Allocate the number of points and polygons. Note, the number of 545 ! polygons is identical to the number of surfaces elements, whereas the 546 ! number of points (vertices), which define the polygons, can be larger. 541 547 ALLOCATE( surfaces%points(3,1:surfaces%npoints) ) 542 548 ALLOCATE( surfaces%polygons(5,1:surfaces%ns) ) 543 549 ! 544 550 ! Note, PARAVIEW expects consecutively ordered points, in order to 545 ! unambiguously identify surfaces. 551 ! unambiguously identify surfaces. 546 552 ! Hence, all PEs should know where they start counting, depending on the 547 ! number of points on the other PE's with lower MPI rank. 553 ! number of points on the other PE's with lower MPI rank. 548 554 #if defined( __parallel ) 549 555 CALL MPI_ALLGATHER( surfaces%npoints, 1, MPI_INTEGER, & … … 554 560 555 561 ! 556 ! After the number of vertices is counted, repeat the loops and define the557 ! vertices. Start with the horizontal default surfaces.562 ! After the number of vertices is counted, repeat the loops and define 563 ! the vertices. Start with the horizontal default surfaces. 558 564 ! First, however, determine the offset where couting of points should be 559 565 ! started, which is the sum of points of all PE's with lower MPI rank. … … 564 570 i = i + 1 565 571 ENDDO 566 572 567 573 surfaces%npoints = 0 568 574 point_index = 1 569 575 npg = 0 570 576 571 577 DO l = 0, 1 572 578 DO m = 1, surf_def_h(0)%ns 573 579 ! 574 ! Determine the indices of the respective grid cell inside the 580 ! Determine the indices of the respective grid cell inside the 575 581 ! topography. 576 582 i = surf_def_h(0)%i(m) + surf_def_h(0)%ioff … … 578 584 k = surf_def_h(0)%k(m) + surf_def_h(0)%koff 579 585 ! 580 ! Check if the vertices that define the surface element are 586 ! Check if the vertices that define the surface element are 581 587 ! already defined, if not, increment the counter. 582 588 IF ( point_index(k,j,i) < 0 ) THEN 583 surfaces%npoints = surfaces%npoints + 1 589 surfaces%npoints = surfaces%npoints + 1 584 590 point_index(k,j,i) = point_index_count 585 point_index_count = point_index_count + 1 591 point_index_count = point_index_count + 1 586 592 surfaces%points(1,surfaces%npoints) = ( i  0.5_wp ) * dx 587 593 surfaces%points(2,surfaces%npoints) = ( j  0.5_wp ) * dy … … 589 595 ENDIF 590 596 IF ( point_index(k,j,i+1) < 0 ) THEN 591 surfaces%npoints = surfaces%npoints + 1 597 surfaces%npoints = surfaces%npoints + 1 592 598 point_index(k,j,i+1) = point_index_count 593 point_index_count = point_index_count + 1 599 point_index_count = point_index_count + 1 594 600 surfaces%points(1,surfaces%npoints) = ( i + 1  0.5_wp ) * dx 595 601 surfaces%points(2,surfaces%npoints) = ( j  0.5_wp ) * dy … … 597 603 ENDIF 598 604 IF ( point_index(k,j+1,i+1) < 0 ) THEN 599 surfaces%npoints = surfaces%npoints + 1 605 surfaces%npoints = surfaces%npoints + 1 600 606 point_index(k,j+1,i+1) = point_index_count 601 point_index_count = point_index_count + 1 607 point_index_count = point_index_count + 1 602 608 surfaces%points(1,surfaces%npoints) = ( i + 1  0.5_wp ) * dx 603 609 surfaces%points(2,surfaces%npoints) = ( j + 1  0.5_wp ) * dy … … 605 611 ENDIF 606 612 IF ( point_index(k,j+1,i) < 0 ) THEN 607 surfaces%npoints = surfaces%npoints + 1 613 surfaces%npoints = surfaces%npoints + 1 608 614 point_index(k,j+1,i) = point_index_count 609 point_index_count = point_index_count + 1 615 point_index_count = point_index_count + 1 610 616 surfaces%points(1,surfaces%npoints) = ( i  0.5_wp ) * dx 611 617 surfaces%points(2,surfaces%npoints) = ( j + 1  0.5_wp ) * dy 612 618 surfaces%points(3,surfaces%npoints) = zw(k) 613 619 ENDIF 614 620 615 621 npg = npg + 1 616 622 surfaces%polygons(1,npg) = 4 … … 626 632 k = surf_lsm_h%k(m) + surf_lsm_h%koff 627 633 IF ( point_index(k,j,i) < 0 ) THEN 628 surfaces%npoints = surfaces%npoints + 1 634 surfaces%npoints = surfaces%npoints + 1 629 635 point_index(k,j,i) = point_index_count 630 point_index_count = point_index_count + 1 636 point_index_count = point_index_count + 1 631 637 surfaces%points(1,surfaces%npoints) = ( i  0.5_wp ) * dx 632 638 surfaces%points(2,surfaces%npoints) = ( j  0.5_wp ) * dy … … 634 640 ENDIF 635 641 IF ( point_index(k,j,i+1) < 0 ) THEN 636 surfaces%npoints = surfaces%npoints + 1 642 surfaces%npoints = surfaces%npoints + 1 637 643 point_index(k,j,i+1) = point_index_count 638 point_index_count = point_index_count + 1 644 point_index_count = point_index_count + 1 639 645 surfaces%points(1,surfaces%npoints) = ( i + 1  0.5_wp ) * dx 640 646 surfaces%points(2,surfaces%npoints) = ( j  0.5_wp ) * dy … … 642 648 ENDIF 643 649 IF ( point_index(k,j+1,i+1) < 0 ) THEN 644 surfaces%npoints = surfaces%npoints + 1 650 surfaces%npoints = surfaces%npoints + 1 645 651 point_index(k,j+1,i+1) = point_index_count 646 point_index_count = point_index_count + 1 652 point_index_count = point_index_count + 1 647 653 surfaces%points(1,surfaces%npoints) = ( i + 1  0.5_wp ) * dx 648 654 surfaces%points(2,surfaces%npoints) = ( j + 1  0.5_wp ) * dy … … 650 656 ENDIF 651 657 IF ( point_index(k,j+1,i) < 0 ) THEN 652 surfaces%npoints = surfaces%npoints + 1 658 surfaces%npoints = surfaces%npoints + 1 653 659 point_index(k,j+1,i) = point_index_count 654 point_index_count = point_index_count + 1 660 point_index_count = point_index_count + 1 655 661 surfaces%points(1,surfaces%npoints) = ( i  0.5_wp ) * dx 656 662 surfaces%points(2,surfaces%npoints) = ( j + 1  0.5_wp ) * dy 657 663 surfaces%points(3,surfaces%npoints) = zw(k) 658 664 ENDIF 659 665 660 666 npg = npg + 1 661 667 surfaces%polygons(1,npg) = 4 … … 663 669 surfaces%polygons(3,npg) = point_index(k,j,i+1) 664 670 surfaces%polygons(4,npg) = point_index(k,j+1,i+1) 665 surfaces%polygons(5,npg) = point_index(k,j+1,i) 671 surfaces%polygons(5,npg) = point_index(k,j+1,i) 666 672 ENDDO 667 673 668 674 DO m = 1, surf_usm_h%ns 669 675 i = surf_usm_h%i(m) + surf_usm_h%ioff 670 676 j = surf_usm_h%j(m) + surf_usm_h%joff 671 677 k = surf_usm_h%k(m) + surf_usm_h%koff 672 678 673 679 IF ( point_index(k,j,i) < 0 ) THEN 674 surfaces%npoints = surfaces%npoints + 1 680 surfaces%npoints = surfaces%npoints + 1 675 681 point_index(k,j,i) = point_index_count 676 point_index_count = point_index_count + 1 682 point_index_count = point_index_count + 1 677 683 surfaces%points(1,surfaces%npoints) = ( i  0.5_wp ) * dx 678 684 surfaces%points(2,surfaces%npoints) = ( j  0.5_wp ) * dy … … 680 686 ENDIF 681 687 IF ( point_index(k,j,i+1) < 0 ) THEN 682 surfaces%npoints = surfaces%npoints + 1 688 surfaces%npoints = surfaces%npoints + 1 683 689 point_index(k,j,i+1) = point_index_count 684 point_index_count = point_index_count + 1 690 point_index_count = point_index_count + 1 685 691 surfaces%points(1,surfaces%npoints) = ( i + 1  0.5_wp ) * dx 686 692 surfaces%points(2,surfaces%npoints) = ( j  0.5_wp ) * dy … … 688 694 ENDIF 689 695 IF ( point_index(k,j+1,i+1) < 0 ) THEN 690 surfaces%npoints = surfaces%npoints + 1 696 surfaces%npoints = surfaces%npoints + 1 691 697 point_index(k,j+1,i+1) = point_index_count 692 point_index_count = point_index_count + 1 698 point_index_count = point_index_count + 1 693 699 surfaces%points(1,surfaces%npoints) = ( i + 1  0.5_wp ) * dx 694 700 surfaces%points(2,surfaces%npoints) = ( j + 1  0.5_wp ) * dy … … 696 702 ENDIF 697 703 IF ( point_index(k,j+1,i) < 0 ) THEN 698 surfaces%npoints = surfaces%npoints + 1 704 surfaces%npoints = surfaces%npoints + 1 699 705 point_index(k,j+1,i) = point_index_count 700 point_index_count = point_index_count + 1 706 point_index_count = point_index_count + 1 701 707 surfaces%points(1,surfaces%npoints) = ( i  0.5_wp ) * dx 702 708 surfaces%points(2,surfaces%npoints) = ( j + 1  0.5_wp ) * dy 703 709 surfaces%points(3,surfaces%npoints) = zw(k) 704 710 ENDIF 705 711 706 712 npg = npg + 1 707 713 surfaces%polygons(1,npg) = 4 … … 709 715 surfaces%polygons(3,npg) = point_index(k,j,i+1) 710 716 surfaces%polygons(4,npg) = point_index(k,j+1,i+1) 711 surfaces%polygons(5,npg) = point_index(k,j+1,i) 717 surfaces%polygons(5,npg) = point_index(k,j+1,i) 712 718 ENDDO 713 719 714 720 DO l = 0, 3 715 721 DO m = 1, surf_def_v(l)%ns 716 ! 717 ! Determine the indices of the respective grid cell inside the topography. 718 ! Please note, jindex for northfacing surfaces ( l==0 ) is 719 ! identical to the reference jindex outside the grid box. 722 ! 723 ! Determine the indices of the respective grid cell inside the 724 ! topography. 725 ! NOTE, jindex for northfacing surfaces ( l==0 ) is 726 ! identical to the reference jindex outside the grid box. 720 727 ! Equivalent for eastfacing surfaces and iindex. 721 728 i = surf_def_v(l)%i(m) + MERGE( surf_def_v(l)%ioff, 0, l == 3 ) 722 729 j = surf_def_v(l)%j(m) + MERGE( surf_def_v(l)%joff, 0, l == 1 ) 723 730 k = surf_def_v(l)%k(m) + surf_def_v(l)%koff 724 ! 731 ! 725 732 ! Lower left /front vertex 726 733 IF ( point_index(k,j,i) < 0 ) THEN 727 surfaces%npoints = surfaces%npoints + 1 734 surfaces%npoints = surfaces%npoints + 1 728 735 point_index(k,j,i) = point_index_count 729 point_index_count = point_index_count + 1 736 point_index_count = point_index_count + 1 730 737 surfaces%points(1,surfaces%npoints) = ( i  0.5_wp ) * dx 731 738 surfaces%points(2,surfaces%npoints) = ( j  0.5_wp ) * dy 732 739 surfaces%points(3,surfaces%npoints) = zw(k1) 733 ENDIF 734 ! 740 ENDIF 741 ! 735 742 ! Upper / lower right index for north and southfacing surfaces 736 743 IF ( l == 0 .OR. l == 1 ) THEN 737 744 IF ( point_index(k,j,i+1) < 0 ) THEN 738 surfaces%npoints = surfaces%npoints + 1 745 surfaces%npoints = surfaces%npoints + 1 739 746 point_index(k,j,i+1) = point_index_count 740 747 point_index_count = point_index_count + 1 … … 742 749 surfaces%points(2,surfaces%npoints) = ( j  0.5_wp ) * dy 743 750 surfaces%points(3,surfaces%npoints) = zw(k1) 744 ENDIF 751 ENDIF 745 752 IF ( point_index(k+1,j,i+1) < 0 ) THEN 746 surfaces%npoints = surfaces%npoints + 1 753 surfaces%npoints = surfaces%npoints + 1 747 754 point_index(k+1,j,i+1) = point_index_count 748 755 point_index_count = point_index_count + 1 … … 751 758 surfaces%points(3,surfaces%npoints) = zw(k) 752 759 ENDIF 753 ! 760 ! 754 761 ! Upper / lower front index for east and westfacing surfaces 755 762 ELSEIF ( l == 2 .OR. l == 3 ) THEN 756 763 IF ( point_index(k,j+1,i) < 0 ) THEN 757 surfaces%npoints = surfaces%npoints + 1 764 surfaces%npoints = surfaces%npoints + 1 758 765 point_index(k,j+1,i) = point_index_count 759 766 point_index_count = point_index_count + 1 … … 761 768 surfaces%points(2,surfaces%npoints) = ( j + 1  0.5_wp ) * dy 762 769 surfaces%points(3,surfaces%npoints) = zw(k1) 763 ENDIF 770 ENDIF 764 771 IF ( point_index(k+1,j+1,i) < 0 ) THEN 765 surfaces%npoints = surfaces%npoints + 1 772 surfaces%npoints = surfaces%npoints + 1 766 773 point_index(k+1,j+1,i) = point_index_count 767 774 point_index_count = point_index_count + 1 … … 771 778 ENDIF 772 779 ENDIF 773 ! 780 ! 774 781 ! Upper left / front vertex 775 782 IF ( point_index(k+1,j,i) < 0 ) THEN 776 surfaces%npoints = surfaces%npoints + 1 783 surfaces%npoints = surfaces%npoints + 1 777 784 point_index(k+1,j,i) = point_index_count 778 785 point_index_count = point_index_count + 1 … … 781 788 surfaces%points(3,surfaces%npoints) = zw(k) 782 789 ENDIF 783 784 npg = npg + 1 790 791 npg = npg + 1 785 792 IF ( l == 0 .OR. l == 1 ) THEN 786 793 surfaces%polygons(1,npg) = 4 … … 788 795 surfaces%polygons(3,npg) = point_index(k,j,i+1) 789 796 surfaces%polygons(4,npg) = point_index(k+1,j,i+1) 790 surfaces%polygons(5,npg) = point_index(k+1,j,i) 797 surfaces%polygons(5,npg) = point_index(k+1,j,i) 791 798 ELSE 792 799 surfaces%polygons(1,npg) = 4 … … 796 803 surfaces%polygons(5,npg) = point_index(k+1,j,i) 797 804 ENDIF 798 805 799 806 ENDDO 800 807 801 808 DO m = 1, surf_lsm_v(l)%ns 802 809 i = surf_lsm_v(l)%i(m) + MERGE( surf_lsm_v(l)%ioff, 0, l == 3 ) 803 810 j = surf_lsm_v(l)%j(m) + MERGE( surf_lsm_v(l)%joff, 0, l == 1 ) 804 811 k = surf_lsm_v(l)%k(m) + surf_lsm_v(l)%koff 805 ! 812 ! 806 813 ! Lower left /front vertex 807 814 IF ( point_index(k,j,i) < 0 ) THEN 808 surfaces%npoints = surfaces%npoints + 1 815 surfaces%npoints = surfaces%npoints + 1 809 816 point_index(k,j,i) = point_index_count 810 817 point_index_count = point_index_count + 1 … … 812 819 surfaces%points(2,surfaces%npoints) = ( j  0.5_wp ) * dy 813 820 surfaces%points(3,surfaces%npoints) = zw(k1) 814 ENDIF 815 ! 821 ENDIF 822 ! 816 823 ! Upper / lower right index for north and southfacing surfaces 817 824 IF ( l == 0 .OR. l == 1 ) THEN 818 825 IF ( point_index(k,j,i+1) < 0 ) THEN 819 surfaces%npoints = surfaces%npoints + 1 826 surfaces%npoints = surfaces%npoints + 1 820 827 point_index(k,j,i+1) = point_index_count 821 828 point_index_count = point_index_count + 1 … … 823 830 surfaces%points(2,surfaces%npoints) = ( j  0.5_wp ) * dy 824 831 surfaces%points(3,surfaces%npoints) = zw(k1) 825 ENDIF 832 ENDIF 826 833 IF ( point_index(k+1,j,i+1) < 0 ) THEN 827 surfaces%npoints = surfaces%npoints + 1 834 surfaces%npoints = surfaces%npoints + 1 828 835 point_index(k+1,j,i+1) = point_index_count 829 836 point_index_count = point_index_count + 1 … … 832 839 surfaces%points(3,surfaces%npoints) = zw(k) 833 840 ENDIF 834 ! 841 ! 835 842 ! Upper / lower front index for east and westfacing surfaces 836 843 ELSEIF ( l == 2 .OR. l == 3 ) THEN 837 844 IF ( point_index(k,j+1,i) < 0 ) THEN 838 surfaces%npoints = surfaces%npoints + 1 845 surfaces%npoints = surfaces%npoints + 1 839 846 point_index(k,j+1,i) = point_index_count 840 847 point_index_count = point_index_count + 1 … … 842 849 surfaces%points(2,surfaces%npoints) = ( j + 1  0.5_wp ) * dy 843 850 surfaces%points(3,surfaces%npoints) = zw(k1) 844 ENDIF 851 ENDIF 845 852 IF ( point_index(k+1,j+1,i) < 0 ) THEN 846 surfaces%npoints = surfaces%npoints + 1 853 surfaces%npoints = surfaces%npoints + 1 847 854 point_index(k+1,j+1,i) = point_index_count 848 855 point_index_count = point_index_count + 1 … … 852 859 ENDIF 853 860 ENDIF 854 ! 861 ! 855 862 ! Upper left / front vertex 856 863 IF ( point_index(k+1,j,i) < 0 ) THEN 857 surfaces%npoints = surfaces%npoints + 1 864 surfaces%npoints = surfaces%npoints + 1 858 865 point_index(k+1,j,i) = point_index_count 859 866 point_index_count = point_index_count + 1 … … 862 869 surfaces%points(3,surfaces%npoints) = zw(k) 863 870 ENDIF 864 865 npg = npg + 1 871 872 npg = npg + 1 866 873 IF ( l == 0 .OR. l == 1 ) THEN 867 874 surfaces%polygons(1,npg) = 4 … … 869 876 surfaces%polygons(3,npg) = point_index(k,j,i+1) 870 877 surfaces%polygons(4,npg) = point_index(k+1,j,i+1) 871 surfaces%polygons(5,npg) = point_index(k+1,j,i) 878 surfaces%polygons(5,npg) = point_index(k+1,j,i) 872 879 ELSE 873 880 surfaces%polygons(1,npg) = 4 … … 882 889 j = surf_usm_v(l)%j(m) + MERGE( surf_usm_v(l)%joff, 0, l == 1 ) 883 890 k = surf_usm_v(l)%k(m) + surf_usm_v(l)%koff 884 ! 891 ! 885 892 ! Lower left /front vertex 886 893 IF ( point_index(k,j,i) < 0 ) THEN 887 surfaces%npoints = surfaces%npoints + 1 894 surfaces%npoints = surfaces%npoints + 1 888 895 point_index(k,j,i) = point_index_count 889 896 point_index_count = point_index_count + 1 … … 891 898 surfaces%points(2,surfaces%npoints) = ( j  0.5_wp ) * dy 892 899 surfaces%points(3,surfaces%npoints) = zw(k1) 893 ENDIF 894 ! 900 ENDIF 901 ! 895 902 ! Upper / lower right index for north and southfacing surfaces 896 903 IF ( l == 0 .OR. l == 1 ) THEN 897 904 IF ( point_index(k,j,i+1) < 0 ) THEN 898 surfaces%npoints = surfaces%npoints + 1 905 surfaces%npoints = surfaces%npoints + 1 899 906 point_index(k,j,i+1) = point_index_count 900 907 point_index_count = point_index_count + 1 … … 902 909 surfaces%points(2,surfaces%npoints) = ( j  0.5_wp ) * dy 903 910 surfaces%points(3,surfaces%npoints) = zw(k1) 904 ENDIF 911 ENDIF 905 912 IF ( point_index(k+1,j,i+1) < 0 ) THEN 906 surfaces%npoints = surfaces%npoints + 1 913 surfaces%npoints = surfaces%npoints + 1 907 914 point_index(k+1,j,i+1) = point_index_count 908 915 point_index_count = point_index_count + 1 … … 911 918 surfaces%points(3,surfaces%npoints) = zw(k) 912 919 ENDIF 913 ! 920 ! 914 921 ! Upper / lower front index for east and westfacing surfaces 915 922 ELSEIF ( l == 2 .OR. l == 3 ) THEN 916 923 IF ( point_index(k,j+1,i) < 0 ) THEN 917 surfaces%npoints = surfaces%npoints + 1 924 surfaces%npoints = surfaces%npoints + 1 918 925 point_index(k,j+1,i) = point_index_count 919 926 point_index_count = point_index_count + 1 … … 921 928 surfaces%points(2,surfaces%npoints) = ( j + 1  0.5_wp ) * dy 922 929 surfaces%points(3,surfaces%npoints) = zw(k1) 923 ENDIF 930 ENDIF 924 931 IF ( point_index(k+1,j+1,i) < 0 ) THEN 925 surfaces%npoints = surfaces%npoints + 1 932 surfaces%npoints = surfaces%npoints + 1 926 933 point_index(k+1,j+1,i) = point_index_count 927 934 point_index_count = point_index_count + 1 … … 931 938 ENDIF 932 939 ENDIF 933 ! 940 ! 934 941 ! Upper left / front vertex 935 942 IF ( point_index(k+1,j,i) < 0 ) THEN 936 surfaces%npoints = surfaces%npoints + 1 943 surfaces%npoints = surfaces%npoints + 1 937 944 point_index(k+1,j,i) = point_index_count 938 945 point_index_count = point_index_count + 1 … … 941 948 surfaces%points(3,surfaces%npoints) = zw(k) 942 949 ENDIF 943 944 npg = npg + 1 950 951 npg = npg + 1 945 952 IF ( l == 0 .OR. l == 1 ) THEN 946 953 surfaces%polygons(1,npg) = 4 … … 948 955 surfaces%polygons(3,npg) = point_index(k,j,i+1) 949 956 surfaces%polygons(4,npg) = point_index(k+1,j,i+1) 950 surfaces%polygons(5,npg) = point_index(k+1,j,i) 957 surfaces%polygons(5,npg) = point_index(k+1,j,i) 951 958 ELSE 952 959 surfaces%polygons(1,npg) = 4 … … 957 964 ENDIF 958 965 ENDDO 959 966 960 967 ENDDO 961 ! 968 ! 962 969 ! Deallocate temporary dummy variable 963 970 DEALLOCATE ( point_index ) 964 971 ! 965 972 ! Sumup total number of vertices on domain. This 966 ! will be needed for postprocessing. 973 ! will be needed for postprocessing. 967 974 surfaces%npoints_total = 0 968 975 #if defined( __parallel ) … … 974 981 ENDIF 975 982 ! 976 ! If output to netcdf is enabled, setup the coordinate arrays that 977 ! unambiguously describe the position and orientation of each surface 978 ! element. 983 ! If output to netcdf is enabled, setup the coordinate arrays that 984 ! unambiguously describe the position and orientation of each surface 985 ! element. 979 986 IF ( to_netcdf ) THEN 980 987 ! … … 990 997 ! 991 998 ! Gather the number of surface on each processor, in order to number 992 ! the surface elements in ascending order with respect to the total 999 ! the surface elements in ascending order with respect to the total 993 1000 ! number of surfaces in the domain. 994 1001 #if defined( __parallel ) … … 999 1006 #endif 1000 1007 ! 1001 ! First, however, determine the offset where couting of the surfaces 1008 ! First, however, determine the offset where couting of the surfaces 1002 1009 ! should start (the sum of surfaces on all PE's with lower MPI rank). 1003 1010 i = 0 … … 1007 1014 i = i + 1 1008 1015 ENDDO 1009 ! 1010 ! Set coordinate arrays. For horizontal surfaces, azimuth 1011 ! angles are not defined (fill value). Zenith angle is 0 (180) for 1016 ! 1017 ! Set coordinate arrays. For horizontal surfaces, azimuth 1018 ! angles are not defined (fill value). Zenith angle is 0 (180) for 1012 1019 ! upward (downward)facing surfaces. 1013 1020 i = start_count … … 1053 1060 mm = mm + 1 1054 1061 ENDDO 1055 ! 1056 ! For vertical surfaces, zenith angles are not defined (fill value). 1057 ! Azimuth angle: eastward (0), westward (180), northward (270), 1058 ! southward (90). 1062 ! 1063 ! For vertical surfaces, zenith angles are not defined (fill value). 1064 ! Azimuth angle: eastward (0), westward (180), northward (270), 1065 ! southward (90). 1059 1066 DO l = 0, 3 1060 1067 IF ( l == 0 ) THEN … … 1075 1082 off_y = 0.0_wp 1076 1083 ENDIF 1077 1084 1078 1085 DO m = 1, surf_def_v(l)%ns 1079 1086 surfaces%s(mm) = i … … 1107 1114 ENDDO 1108 1115 ENDDO 1109 ! 1110 ! Finally, define UTM coordinates, which are the x/ycoordinates 1116 ! 1117 ! Finally, define UTM coordinates, which are the x/ycoordinates 1111 1118 ! plus the origin (lowerleft coordinate of the model domain). 1112 surfaces%es_utm = surfaces%xs + init_model%origin_x 1113 surfaces%ns_utm = surfaces%ys + init_model%origin_y 1114 ! 1115 ! Initialize NetCDF data output. Please note, local start position for 1116 ! the surface elements in the NetCDF file is surfaces%s(1), while 1119 surfaces%es_utm = surfaces%xs + init_model%origin_x 1120 surfaces%ns_utm = surfaces%ys + init_model%origin_y 1121 ! 1122 ! Initialize NetCDF data output. Please note, local start position for 1123 ! the surface elements in the NetCDF file is surfaces%s(1), while 1117 1124 ! the number of surfaces on the subdomain is given by surfaces%ns. 1118 1125 #if defined( __netcdf4_parallel ) … … 1120 1127 ! 1121 1128 ! Calculate number of time steps to be output 1122 ntdim_surf(0) = dosurf_time_count(0) + CEILING( & 1123 ( end_time  MAX( & 1124 MERGE(skip_time_dosurf, skip_time_dosurf + spinup_time, & 1125 data_output_during_spinup ), & 1126 simulated_time_at_begin ) & 1129 ntdim_surf(0) = dosurf_time_count(0) + CEILING( & 1130 ( end_time  MAX( & 1131 MERGE( skip_time_dosurf, & 1132 skip_time_dosurf + spinup_time, & 1133 data_output_during_spinup ), & 1134 simulated_time_at_begin ) & 1127 1135 ) / dt_dosurf ) 1128 1136 1129 1137 ntdim_surf(1) = dosurf_time_count(1) + CEILING( & 1130 1138 ( end_time  MAX( & 1131 MERGE( skip_time_dosurf_av, skip_time_dosurf_av & 1132 + spinup_time, data_output_during_spinup ), & 1139 MERGE( skip_time_dosurf_av, & 1140 skip_time_dosurf_av + spinup_time, & 1141 data_output_during_spinup ), & 1133 1142 simulated_time_at_begin ) & 1134 1143 ) / dt_dosurf_av ) … … 1152 1161 IOR( NF90_NOCLOBBER, & 1153 1162 IOR( NF90_NETCDF4, NF90_MPIIO ) ), & 1154 id_set_surf(av), COMM = comm2d, INFO = MPI_INFO_NULL ) 1163 id_set_surf(av), & 1164 COMM = comm2d, INFO = MPI_INFO_NULL ) 1155 1165 CALL netcdf_handle_error( 'surface_data_output_mod', 5550 ) 1156 1166 1157 1167 ! Write some global attributes 1158 1168 IF ( av == 0 ) THEN 1159 CALL netcdf_create_global_atts( id_set_surf(av), 'surfacedata', TRIM( run_description_header ), 5551 ) 1169 CALL netcdf_create_global_atts( id_set_surf(av), & 1170 'surfacedata', & 1171 TRIM( run_description_header ),& 1172 5551 ) 1160 1173 time_average_text = ' ' 1161 1174 ELSE 1162 CALL netcdf_create_global_atts( id_set_surf(av), 'surfacedata_av', TRIM( run_description_header ), 5552 ) 1163 WRITE ( time_average_text,'(F7.1,'' s avg'')' ) averaging_interval_surf 1164 nc_stat = NF90_PUT_ATT( id_set_surf(av), NF90_GLOBAL, 'time_avg', & 1175 CALL netcdf_create_global_atts( id_set_surf(av), & 1176 'surfacedata_av', & 1177 TRIM( run_description_header ),& 1178 5552 ) 1179 WRITE ( time_average_text,'(F7.1,'' s avg'')' ) & 1180 averaging_interval_surf 1181 nc_stat = NF90_PUT_ATT( id_set_surf(av), NF90_GLOBAL, & 1182 'time_avg', & 1165 1183 TRIM( time_average_text ) ) 1166 1184 CALL netcdf_handle_error( 'surface_data_output_mod', 5553 ) … … 1170 1188 ! 1171 1189 ! Define time coordinate for surface data. 1172 ! For parallel output the time dimension has to be limited (ntdim_surf),1173 ! otherwise the performance drops significantly.1190 ! For parallel output the time dimension has to be limited 1191 ! (ntdim_surf), otherwise the performance drops significantly. 1174 1192 CALL netcdf_create_dim( id_set_surf(av), 'time', ntdim_surf(av), & 1175 1193 id_dim_time_surf(av), 5554 ) 1176 1194 1177 CALL netcdf_create_var( id_set_surf(av), (/ id_dim_time_surf(av) /), & 1178 'time', NF90_DOUBLE, id_var_time_surf(av), & 1179 'seconds since '//TRIM(init_model%origin_time), & 1195 CALL netcdf_create_var( id_set_surf(av), & 1196 (/ id_dim_time_surf(av) /), & 1197 'time', NF90_DOUBLE, & 1198 id_var_time_surf(av), & 1199 'seconds since '// & 1200 TRIM(init_model%origin_time), & 1180 1201 'time', 5555, 5555, 5555 ) 1181 CALL netcdf_create_att( id_set_surf(av), id_var_time_surf(av), 'standard_name', 'time', 5556) 1182 CALL netcdf_create_att( id_set_surf(av), id_var_time_surf(av), 'axis', 'T', 5557) 1202 1203 CALL netcdf_create_att( id_set_surf(av), id_var_time_surf(av), & 1204 'standard_name', 'time', 5556) 1205 1206 CALL netcdf_create_att( id_set_surf(av), id_var_time_surf(av), & 1207 'axis', 'T', 5557) 1183 1208 ! 1184 1209 ! Define spatial dimensions and coordinates: … … 1188 1213 CALL netcdf_create_var( id_set_surf(av), (/ id_dim_s_surf(av) /), & 1189 1214 's', NF90_DOUBLE, id_var_s_surf(av), & 1190 '1', '', 5559, 5559, 5559 ) 1215 '1', 'number of surface element', & 1216 5559, 5559, 5559 ) 1191 1217 ! 1192 1218 ! Define x coordinate 1193 1219 CALL netcdf_create_var( id_set_surf(av), (/ id_dim_s_surf(av) /), & 1194 1220 'xs', NF90_DOUBLE, id_var_xs_surf(av), & 1195 'meters', '', 5561, 5561, 5561 ) 1221 'meters', & 1222 'distance to origin in xdirection', & 1223 5561, 5561, 5561 ) 1196 1224 ! 1197 1225 ! Define y coordinate 1198 1226 CALL netcdf_create_var( id_set_surf(av), (/ id_dim_s_surf(av) /), & 1199 1227 'ys', NF90_DOUBLE, id_var_ys_surf(av), & 1200 'meters', '', 5562, 5562, 5562 ) 1228 'meters', & 1229 'distance to origin in ydirection', & 1230 5562, 5562, 5562 ) 1201 1231 ! 1202 1232 ! Define z coordinate 1203 1233 CALL netcdf_create_var( id_set_surf(av), (/ id_dim_s_surf(av) /), & 1204 1234 'zs', NF90_DOUBLE, id_var_zs_surf(av), & 1205 'meters', '', 5560, 5560, 5560 ) 1235 'meters', 'height', 5560, 5560, 5560 ) 1236 CALL netcdf_create_att( id_set_surf(av), id_var_zs_surf(av), & 1237 'standard_name', 'height', 5583 ) 1206 1238 1207 1239 ! 1208 1240 ! Define UTM coordinates 1209 CALL netcdf_create_var( id_set_surf(av), (/ id_dim_s_surf(av) /), & 1210 'Es_UTM', NF90_DOUBLE, id_var_etum_surf(av), & 1241 CALL netcdf_create_var( id_set_surf(av), (/ id_dim_s_surf(av) /), & 1242 'Es_UTM', NF90_DOUBLE, & 1243 id_var_etum_surf(av), & 1211 1244 'meters', '', 5563, 5563, 5563 ) 1212 CALL netcdf_create_var( id_set_surf(av), (/ id_dim_s_surf(av) /), & 1213 'Ns_UTM', NF90_DOUBLE, id_var_nutm_surf(av), & 1245 CALL netcdf_create_var( id_set_surf(av), (/ id_dim_s_surf(av) /), & 1246 'Ns_UTM', NF90_DOUBLE, & 1247 id_var_nutm_surf(av), & 1214 1248 'meters', '', 5564, 5564, 5564 ) 1249 1250 ! 1251 ! Define angles 1252 CALL netcdf_create_var( id_set_surf(av), (/ id_dim_s_surf(av) /), & 1253 'azimuth', NF90_DOUBLE, & 1254 id_var_azimuth_surf(av), & 1255 'degree', 'azimuth angle', & 1256 5577, 5578, 5579, & 1257 fill = .TRUE. ) 1258 CALL netcdf_create_att( id_set_surf(av), id_var_azimuth_surf(av), & 1259 'standard_name', 'surface_azimuth_angle', & 1260 5584 ) 1261 1262 CALL netcdf_create_var( id_set_surf(av), (/ id_dim_s_surf(av) /), & 1263 'zenith', NF90_DOUBLE, & 1264 id_var_zenith_surf(av), & 1265 'degree', '', 5580, 5581, 5582, & 1266 fill = .TRUE. ) 1215 1267 ! 1216 1268 ! Define the variables … … 1227 1279 ! 1228 1280 ! Set no fill for every variable to increase performance. 1229 nc_stat = NF90_DEF_VAR_FILL( id_set_surf(av), &1230 id_var_dosurf(av,i), &1281 nc_stat = NF90_DEF_VAR_FILL( id_set_surf(av), & 1282 id_var_dosurf(av,i), & 1231 1283 1, 0 ) 1232 1284 CALL netcdf_handle_error( 'surface_data_output_init', 5566 ) 1233 1285 ! 1234 1286 ! Set collective io operations for parallel io 1235 nc_stat = NF90_VAR_PAR_ACCESS( id_set_surf(av), &1236 id_var_dosurf(av,i), &1287 nc_stat = NF90_VAR_PAR_ACCESS( id_set_surf(av), & 1288 id_var_dosurf(av,i), & 1237 1289 NF90_COLLECTIVE ) 1238 1290 CALL netcdf_handle_error( 'surface_data_output_init', 5567 ) … … 1250 1302 1251 1303 ! 1252 ! Set general no fill, otherwise the performance drops significantly for1253 ! parallel output.1304 ! Set general no fill, otherwise the performance drops significantly 1305 ! for parallel output. 1254 1306 nc_stat = NF90_SET_FILL( id_set_surf(av), NF90_NOFILL, oldmode ) 1255 1307 CALL netcdf_handle_error( 'surface_data_output_init', 5569 ) … … 1272 1324 ENDDO 1273 1325 1274 nc_stat = NF90_PUT_VAR( id_set_surf(av), id_var_s_surf(av), &1275 netcdf_data_1d, start = (/ 1 /), &1326 nc_stat = NF90_PUT_VAR( id_set_surf(av), id_var_s_surf(av), & 1327 netcdf_data_1d, start = (/ 1 /), & 1276 1328 count = (/ surfaces%ns_total /) ) 1277 1329 CALL netcdf_handle_error( 'surface_data_output_init', 5571 ) … … 1298 1350 CALL netcdf_handle_error( 'surface_data_output_init', 5574 ) 1299 1351 1300 nc_stat = NF90_PUT_VAR( id_set_surf(av), id_var_etum_surf(av), & 1301 surfaces%es_utm, start = (/ surfaces%s(1) /), & 1352 nc_stat = NF90_PUT_VAR( id_set_surf(av), id_var_etum_surf(av), & 1353 surfaces%es_utm, & 1354 start = (/ surfaces%s(1) /), & 1302 1355 count = (/ surfaces%ns /) ) 1303 1356 CALL netcdf_handle_error( 'surface_data_output_init', 5575 ) 1304 1357 1305 nc_stat = NF90_PUT_VAR( id_set_surf(av), id_var_nutm_surf(av), & 1306 surfaces%ns_utm, start = (/ surfaces%s(1) /), & 1358 nc_stat = NF90_PUT_VAR( id_set_surf(av), id_var_nutm_surf(av), & 1359 surfaces%ns_utm, & 1360 start = (/ surfaces%s(1) /), & 1307 1361 count = (/ surfaces%ns /) ) 1308 1362 CALL netcdf_handle_error( 'surface_data_output_init', 5576 ) 1309 1363 1364 nc_stat = NF90_PUT_VAR( id_set_surf(av), id_var_azimuth_surf(av), & 1365 surfaces%azimuth, & 1366 start = (/ surfaces%s(1) /), & 1367 count = (/ surfaces%ns /) ) 1368 CALL netcdf_handle_error( 'surface_data_output_init', 5585 ) 1369 1370 nc_stat = NF90_PUT_VAR( id_set_surf(av), id_var_zenith_surf(av), & 1371 surfaces%zenith, & 1372 start = (/ surfaces%s(1) /), & 1373 count = (/ surfaces%ns /) ) 1374 CALL netcdf_handle_error( 'surface_data_output_init', 5586 ) 1375 1310 1376 ENDDO 1311 1377 #endif … … 1314 1380 1315 1381 END SUBROUTINE surface_data_output_init 1316 1382 1317 1383 !! 1318 1384 ! Description: 1319 1385 !  1320 !> Routine for controlling the data output. Surface data is collected from 1321 !> different types of surfaces (default, natural, urban) and different 1386 !> Routine for controlling the data output. Surface data is collected from 1387 !> different types of surfaces (default, natural, urban) and different 1322 1388 !> orientation and written to one 1Doutput array. Further, NetCDF routines 1323 1389 !> are called to write the surface data in the respective NetCDF files. 1324 !! 1390 !! 1325 1391 SUBROUTINE surface_data_output( av ) 1326 1392 1327 1393 USE control_parameters, & 1328 1394 ONLY: io_blocks, io_group, time_since_reference_point … … 1331 1397 ONLY: comm2d, ierr 1332 1398 1333 1399 1334 1400 IMPLICIT NONE 1335 1401 1336 1402 CHARACTER(LEN=100) :: trimvar = ' ' !< dummy for single output variable 1337 1403 1338 1404 INTEGER(iwp) :: av !< id indicating average or nonaverage data output 1339 1405 INTEGER(iwp) :: i !< loop index … … 1344 1410 IF ( dosurf_no(av) == 0 ) RETURN 1345 1411 ! 1346 ! In case of VTK output, check if binary files are open and write coordinates. 1412 ! In case of VTK output, check if binary files are open and write coordinates. 1347 1413 IF ( to_vtk ) THEN 1348 1414 … … 1372 1438 #if defined( __netcdf4_parallel ) 1373 1439 IF ( dosurf_time_count(av) + 1 > ntdim_surf(av) ) THEN 1374 WRITE ( message_string, * ) 'Output of surface data is not given at t=', & 1375 time_since_reference_point, 's because the maximum ', & 1376 'number of output time levels is exceeded.' 1440 WRITE ( message_string, * ) & 1441 'Output of surface data is not given at t=', & 1442 time_since_reference_point, 's because the maximum ', & 1443 'number of output time levels is exceeded.' 1377 1444 CALL message( 'surface_data_output', 'PA0539', 0, 1, 0, 6, 0 ) 1378 1445 1379 1446 RETURN 1380 1447 … … 1403 1470 trimvar = TRIM( dosurf(av,n_out) ) 1404 1471 ! 1405 ! Set the output array to the _FillValue in case it is not 1472 ! Set the output array to the _FillValue in case it is not 1406 1473 ! defined for each type of surface. 1407 1474 surfaces%var_out = surfaces%fillvalue … … 1427 1494 surf_def_v(3)%us, & 1428 1495 surf_lsm_v(3)%us, & 1429 surf_usm_v(3)%us ) 1430 ELSE 1431 ! 1432 ! Output of averaged data 1433 surfaces%var_out(:) = surfaces%var_av(:,n_out) / & 1434 REAL( average_count_surf, KIND=wp ) 1435 surfaces%var_av(:,n_out) = 0.0_wp 1436 1496 surf_usm_v(3)%us ) 1497 ELSE 1498 ! 1499 ! Output of averaged data 1500 surfaces%var_out(:) = surfaces%var_av(:,n_out) / & 1501 REAL( average_count_surf, KIND=wp ) 1502 surfaces%var_av(:,n_out) = 0.0_wp 1503 1437 1504 ENDIF 1438 1505 … … 1456 1523 surf_def_v(3)%ts, & 1457 1524 surf_lsm_v(3)%ts, & 1458 surf_usm_v(3)%ts ) 1459 ELSE 1460 ! 1461 ! Output of averaged data 1462 surfaces%var_out(:) = surfaces%var_av(:,n_out) / & 1463 REAL( average_count_surf, KIND=wp ) 1464 surfaces%var_av(:,n_out) = 0.0_wp 1465 1525 surf_usm_v(3)%ts ) 1526 ELSE 1527 ! 1528 ! Output of averaged data 1529 surfaces%var_out(:) = surfaces%var_av(:,n_out) / & 1530 REAL( average_count_surf, KIND=wp ) 1531 surfaces%var_av(:,n_out) = 0.0_wp 1532 1466 1533 ENDIF 1467 1534 … … 1485 1552 surf_def_v(3)%qs, & 1486 1553 surf_lsm_v(3)%qs, & 1487 surf_usm_v(3)%qs ) 1488 ELSE 1489 ! 1490 ! Output of averaged data 1491 surfaces%var_out(:) = surfaces%var_av(:,n_out) / & 1492 REAL( average_count_surf, KIND=wp ) 1493 surfaces%var_av(:,n_out) = 0.0_wp 1494 1554 surf_usm_v(3)%qs ) 1555 ELSE 1556 ! 1557 ! Output of averaged data 1558 surfaces%var_out(:) = surfaces%var_av(:,n_out) / & 1559 REAL( average_count_surf, KIND=wp ) 1560 surfaces%var_av(:,n_out) = 0.0_wp 1561 1495 1562 ENDIF 1496 1563 … … 1514 1581 surf_def_v(3)%ss, & 1515 1582 surf_lsm_v(3)%ss, & 1516 surf_usm_v(3)%ss ) 1517 ELSE 1518 ! 1519 ! Output of averaged data 1520 surfaces%var_out(:) = surfaces%var_av(:,n_out) / & 1521 REAL( average_count_surf, KIND=wp ) 1522 surfaces%var_av(:,n_out) = 0.0_wp 1523 1583 surf_usm_v(3)%ss ) 1584 ELSE 1585 ! 1586 ! Output of averaged data 1587 surfaces%var_out(:) = surfaces%var_av(:,n_out) / & 1588 REAL( average_count_surf, KIND=wp ) 1589 surfaces%var_av(:,n_out) = 0.0_wp 1590 1524 1591 ENDIF 1525 1592 … … 1543 1610 surf_def_v(3)%qcs, & 1544 1611 surf_lsm_v(3)%qcs, & 1545 surf_usm_v(3)%qcs ) 1546 ELSE 1547 ! 1548 ! Output of averaged data 1549 surfaces%var_out(:) = surfaces%var_av(:,n_out) / & 1550 REAL( average_count_surf, KIND=wp ) 1551 surfaces%var_av(:,n_out) = 0.0_wp 1552 1612 surf_usm_v(3)%qcs ) 1613 ELSE 1614 ! 1615 ! Output of averaged data 1616 surfaces%var_out(:) = surfaces%var_av(:,n_out) / & 1617 REAL( average_count_surf, KIND=wp ) 1618 surfaces%var_av(:,n_out) = 0.0_wp 1619 1553 1620 ENDIF 1554 1621 … … 1572 1639 surf_def_v(3)%ncs, & 1573 1640 surf_lsm_v(3)%ncs, & 1574 surf_usm_v(3)%ncs ) 1575 ELSE 1576 ! 1577 ! Output of averaged data 1578 surfaces%var_out(:) = surfaces%var_av(:,n_out) / & 1579 REAL( average_count_surf, KIND=wp ) 1580 surfaces%var_av(:,n_out) = 0.0_wp 1581 1641 surf_usm_v(3)%ncs ) 1642 ELSE 1643 ! 1644 ! Output of averaged data 1645 surfaces%var_out(:) = surfaces%var_av(:,n_out) / & 1646 REAL( average_count_surf, KIND=wp ) 1647 surfaces%var_av(:,n_out) = 0.0_wp 1648 1582 1649 ENDIF 1583 1650 … … 1601 1668 surf_def_v(3)%qrs, & 1602 1669 surf_lsm_v(3)%qrs, & 1603 surf_usm_v(3)%qrs ) 1604 ELSE 1605 ! 1606 ! Output of averaged data 1607 surfaces%var_out(:) = surfaces%var_av(:,n_out) / & 1608 REAL( average_count_surf, KIND=wp ) 1609 surfaces%var_av(:,n_out) = 0.0_wp 1610 1670 surf_usm_v(3)%qrs ) 1671 ELSE 1672 ! 1673 ! Output of averaged data 1674 surfaces%var_out(:) = surfaces%var_av(:,n_out) / & 1675 REAL( average_count_surf, KIND=wp ) 1676 surfaces%var_av(:,n_out) = 0.0_wp 1677 1611 1678 ENDIF 1612 1679 … … 1630 1697 surf_def_v(3)%nrs, & 1631 1698 surf_lsm_v(3)%nrs, & 1632 surf_usm_v(3)%nrs ) 1633 ELSE 1634 ! 1635 ! Output of averaged data 1636 surfaces%var_out(:) = surfaces%var_av(:,n_out) / & 1637 REAL( average_count_surf, KIND=wp ) 1638 surfaces%var_av(:,n_out) = 0.0_wp 1639 1699 surf_usm_v(3)%nrs ) 1700 ELSE 1701 ! 1702 ! Output of averaged data 1703 surfaces%var_out(:) = surfaces%var_av(:,n_out) / & 1704 REAL( average_count_surf, KIND=wp ) 1705 surfaces%var_av(:,n_out) = 0.0_wp 1706 1640 1707 ENDIF 1641 1708 … … 1659 1726 surf_def_v(3)%ol, & 1660 1727 surf_lsm_v(3)%ol, & 1661 surf_usm_v(3)%ol ) 1662 ELSE 1663 ! 1664 ! Output of averaged data 1665 surfaces%var_out(:) = surfaces%var_av(:,n_out) / & 1666 REAL( average_count_surf, KIND=wp ) 1667 surfaces%var_av(:,n_out) = 0.0_wp 1668 1728 surf_usm_v(3)%ol ) 1729 ELSE 1730 ! 1731 ! Output of averaged data 1732 surfaces%var_out(:) = surfaces%var_av(:,n_out) / & 1733 REAL( average_count_surf, KIND=wp ) 1734 surfaces%var_av(:,n_out) = 0.0_wp 1735 1669 1736 ENDIF 1670 1737 … … 1695 1762 REAL( average_count_surf, KIND=wp ) 1696 1763 surfaces%var_av(:,n_out) = 0.0_wp 1697 1764 1698 1765 ENDIF 1699 1766 … … 1724 1791 REAL( average_count_surf, KIND=wp ) 1725 1792 surfaces%var_av(:,n_out) = 0.0_wp 1726 1793 1727 1794 ENDIF 1728 1795 … … 1746 1813 surf_def_v(3)%z0q, & 1747 1814 surf_lsm_v(3)%z0q, & 1748 surf_usm_v(3)%z0q ) 1815 surf_usm_v(3)%z0q ) 1749 1816 ELSE 1750 1817 ! … … 1782 1849 REAL( average_count_surf, KIND=wp ) 1783 1850 surfaces%var_av(:,n_out) = 0.0_wp 1784 1851 1785 1852 ENDIF 1786 1853 … … 1811 1878 REAL( average_count_surf, KIND=wp ) 1812 1879 surfaces%var_av(:,n_out) = 0.0_wp 1813 1880 1814 1881 ENDIF 1815 1882 … … 1840 1907 REAL( average_count_surf, KIND=wp ) 1841 1908 surfaces%var_av(:,n_out) = 0.0_wp 1842 1909 1843 1910 ENDIF 1844 1911 … … 1869 1936 REAL( average_count_surf, KIND=wp ) 1870 1937 surfaces%var_av(:,n_out) = 0.0_wp 1871 1938 1872 1939 ENDIF 1873 1940 … … 1898 1965 REAL( average_count_surf, KIND=wp ) 1899 1966 surfaces%var_av(:,n_out) = 0.0_wp 1900 1967 1901 1968 ENDIF 1902 1969 … … 1955 2022 REAL( average_count_surf, KIND=wp ) 1956 2023 surfaces%var_av(:,n_out) = 0.0_wp 1957 2024 1958 2025 ENDIF 1959 2026 … … 1984 2051 REAL( average_count_surf, KIND=wp ) 1985 2052 surfaces%var_av(:,n_out) = 0.0_wp 1986 2053 1987 2054 ENDIF 1988 2055 … … 2013 2080 REAL( average_count_surf, KIND=wp ) 2014 2081 surfaces%var_av(:,n_out) = 0.0_wp 2015 2082 2016 2083 ENDIF 2017 2084 … … 2042 2109 REAL( average_count_surf, KIND=wp ) 2043 2110 surfaces%var_av(:,n_out) = 0.0_wp 2044 2111 2045 2112 ENDIF 2046 2113 … … 2071 2138 REAL( average_count_surf, KIND=wp ) 2072 2139 surfaces%var_av(:,n_out) = 0.0_wp 2073 2140 2074 2141 ENDIF 2075 2142 … … 2100 2167 REAL( average_count_surf, KIND=wp ) 2101 2168 surfaces%var_av(:,n_out) = 0.0_wp 2102 2169 2103 2170 ENDIF 2104 2171 … … 2129 2196 REAL( average_count_surf, KIND=wp ) 2130 2197 surfaces%var_av(:,n_out) = 0.0_wp 2131 2198 2132 2199 ENDIF 2133 2200 … … 2158 2225 REAL( average_count_surf, KIND=wp ) 2159 2226 surfaces%var_av(:,n_out) = 0.0_wp 2160 2227 2161 2228 ENDIF 2162 2229 … … 2187 2254 REAL( average_count_surf, KIND=wp ) 2188 2255 surfaces%var_av(:,n_out) = 0.0_wp 2189 2256 2190 2257 ENDIF 2191 2258 … … 2216 2283 REAL( average_count_surf, KIND=wp ) 2217 2284 surfaces%var_av(:,n_out) = 0.0_wp 2218 2285 2219 2286 ENDIF 2220 2287 … … 2238 2305 surf_def_v(3)%rad_net, & 2239 2306 surf_lsm_v(3)%rad_net, & 2240 surf_usm_v(3)%rad_net ) 2241 ELSE 2242 ! 2243 ! Output of averaged data 2244 surfaces%var_out(:) = surfaces%var_av(:,n_out) / & 2245 REAL( average_count_surf, KIND=wp ) 2246 surfaces%var_av(:,n_out) = 0.0_wp 2247 2307 surf_usm_v(3)%rad_net ) 2308 ELSE 2309 ! 2310 ! Output of averaged data 2311 surfaces%var_out(:) = surfaces%var_av(:,n_out) / & 2312 REAL( average_count_surf, KIND=wp ) 2313 surfaces%var_av(:,n_out) = 0.0_wp 2314 2248 2315 ENDIF 2249 2316 … … 2274 2341 REAL( average_count_surf, KIND=wp ) 2275 2342 surfaces%var_av(:,n_out) = 0.0_wp 2276 2343 2277 2344 ENDIF 2278 2345 … … 2303 2370 REAL( average_count_surf, KIND=wp ) 2304 2371 surfaces%var_av(:,n_out) = 0.0_wp 2305 2372 2306 2373 ENDIF 2307 2374 … … 2332 2399 REAL( average_count_surf, KIND=wp ) 2333 2400 surfaces%var_av(:,n_out) = 0.0_wp 2334 2401 2335 2402 ENDIF 2336 2403 … … 2361 2428 REAL( average_count_surf, KIND=wp ) 2362 2429 surfaces%var_av(:,n_out) = 0.0_wp 2363 2430 2364 2431 ENDIF 2365 2432 … … 2383 2450 surf_def_v(3)%ghf, & 2384 2451 surf_lsm_v(3)%ghf, & 2385 surf_usm_v(3)%ghf ) 2452 surf_usm_v(3)%ghf ) 2386 2453 ELSE 2387 2454 ! … … 2390 2457 REAL( average_count_surf, KIND=wp ) 2391 2458 surfaces%var_av(:,n_out) = 0.0_wp 2392 2393 ENDIF 2394 2395 CASE ( 'r_a' ) 2459 2460 ENDIF 2461 2462 CASE ( 'r_a' ) 2396 2463 ! 2397 2464 ! Output of instantaneous data … … 2412 2479 surf_def_v(3)%r_a, & 2413 2480 surf_lsm_v(3)%r_a, & 2414 surf_usm_v(3)%r_a ) 2481 surf_usm_v(3)%r_a ) 2415 2482 ELSE 2416 2483 ! … … 2419 2486 REAL( average_count_surf, KIND=wp ) 2420 2487 surfaces%var_av(:,n_out) = 0.0_wp 2421 2422 ENDIF 2423 2424 CASE ( 'r_soil' ) 2488 2489 ENDIF 2490 2491 CASE ( 'r_soil' ) 2425 2492 ! 2426 2493 ! Output of instantaneous data … … 2441 2508 surf_def_v(3)%r_soil, & 2442 2509 surf_lsm_v(3)%r_soil, & 2443 surf_usm_v(3)%r_soil ) 2444 ELSE 2445 ! 2446 ! Output of averaged data 2447 surfaces%var_out(:) = surfaces%var_av(:,n_out) / & 2448 REAL( average_count_surf, KIND=wp ) 2449 surfaces%var_av(:,n_out) = 0.0_wp 2450 2510 surf_usm_v(3)%r_soil ) 2511 ELSE 2512 ! 2513 ! Output of averaged data 2514 surfaces%var_out(:) = surfaces%var_av(:,n_out) / & 2515 REAL( average_count_surf, KIND=wp ) 2516 surfaces%var_av(:,n_out) = 0.0_wp 2517 2451 2518 ENDIF 2452 2519 … … 2470 2537 surf_def_v(3)%r_canopy, & 2471 2538 surf_lsm_v(3)%r_canopy, & 2472 surf_usm_v(3)%r_canopy ) 2473 ELSE 2474 ! 2475 ! Output of averaged data 2476 surfaces%var_out(:) = surfaces%var_av(:,n_out) / & 2477 REAL( average_count_surf, KIND=wp ) 2478 surfaces%var_av(:,n_out) = 0.0_wp 2479 2539 surf_usm_v(3)%r_canopy ) 2540 ELSE 2541 ! 2542 ! Output of averaged data 2543 surfaces%var_out(:) = surfaces%var_av(:,n_out) / & 2544 REAL( average_count_surf, KIND=wp ) 2545 surfaces%var_av(:,n_out) = 0.0_wp 2546 2480 2547 ENDIF 2481 2548 … … 2499 2566 surf_def_v(3)%r_s, & 2500 2567 surf_lsm_v(3)%r_s, & 2501 surf_usm_v(3)%r_s ) 2502 ELSE 2503 ! 2504 ! Output of averaged data 2505 surfaces%var_out(:) = surfaces%var_av(:,n_out) / & 2506 REAL( average_count_surf, KIND=wp ) 2507 surfaces%var_av(:,n_out) = 0.0_wp 2508 2568 surf_usm_v(3)%r_s ) 2569 ELSE 2570 ! 2571 ! Output of averaged data 2572 surfaces%var_out(:) = surfaces%var_av(:,n_out) / & 2573 REAL( average_count_surf, KIND=wp ) 2574 surfaces%var_av(:,n_out) = 0.0_wp 2575 2509 2576 ENDIF 2510 2577 … … 2535 2602 REAL( average_count_surf, KIND=wp ) 2536 2603 surfaces%var_av(:,n_out) = 0.0_wp 2537 2604 2538 2605 ENDIF 2539 2606 … … 2564 2631 REAL( average_count_surf, KIND=wp ) 2565 2632 surfaces%var_av(:,n_out) = 0.0_wp 2566 2633 2567 2634 ENDIF 2568 2635 … … 2593 2660 REAL( average_count_surf, KIND=wp ) 2594 2661 surfaces%var_av(:,n_out) = 0.0_wp 2595 2662 2596 2663 ENDIF 2597 2664 … … 2622 2689 REAL( average_count_surf, KIND=wp ) 2623 2690 surfaces%var_av(:,n_out) = 0.0_wp 2624 2691 2625 2692 ENDIF 2626 2693 … … 2651 2718 REAL( average_count_surf, KIND=wp ) 2652 2719 surfaces%var_av(:,n_out) = 0.0_wp 2653 2720 2654 2721 ENDIF 2655 2722 … … 2680 2747 REAL( average_count_surf, KIND=wp ) 2681 2748 surfaces%var_av(:,n_out) = 0.0_wp 2682 2749 2683 2750 ENDIF 2684 2751 … … 2709 2776 REAL( average_count_surf, KIND=wp ) 2710 2777 surfaces%var_av(:,n_out) = 0.0_wp 2711 2712 ENDIF 2713 2778 2779 ENDIF 2780 2714 2781 CASE ( 'uvw1' ) 2715 2782 ! … … 2738 2805 REAL( average_count_surf, KIND=wp ) 2739 2806 surfaces%var_av(:,n_out) = 0.0_wp 2740 2741 ENDIF 2807 2808 ENDIF 2742 2809 2743 2810 ! … … 2747 2814 END SELECT 2748 2815 ! 2749 ! Write to binary file: 2816 ! Write to binary file: 2750 2817 !  surfaces%points ( 3, 1npoints ) 2751 2818 !  surfaces%polygons ( 5, 1ns ) … … 2768 2835 ENDDO 2769 2836 ENDIF 2770 2837 2771 2838 IF ( to_netcdf ) THEN 2772 2839 #if defined( __netcdf4_parallel ) 2773 2840 ! 2774 2841 ! Write output array to file 2775 nc_stat = NF90_PUT_VAR( id_set_surf(av), id_var_dosurf(av,n_out), & 2776 surfaces%var_out, & 2777 start = (/ surfaces%s(1), dosurf_time_count(av) /), & 2842 nc_stat = NF90_PUT_VAR( id_set_surf(av), id_var_dosurf(av,n_out), & 2843 surfaces%var_out, & 2844 start = (/ surfaces%s(1), & 2845 dosurf_time_count(av) /), & 2778 2846 count = (/ surfaces%ns, 1 /) ) 2779 2847 CALL netcdf_handle_error( 'surface_data_output', 6667 ) … … 2782 2850 2783 2851 ENDDO 2784 2852 2785 2853 ! 2786 2854 ! If averaged output was written to NetCDF file, set the counter to zero 2787 2855 IF ( av == 1 ) average_count_surf = 0 2788 2856 2789 2857 END SUBROUTINE surface_data_output 2790 2858 2791 2859 !! 2792 2860 ! Description: 2793 2861 !  2794 !> Routine for controlling the data averaging. 2795 !! 2862 !> Routine for controlling the data averaging. 2863 !! 2796 2864 SUBROUTINE surface_data_output_averaging 2797 2865 2798 2866 IMPLICIT NONE 2799 2867 2800 2868 CHARACTER(LEN=100) :: trimvar !< dummy variable for current output variable 2801 2869 2802 2870 INTEGER(iwp) :: n_out !< counter variables for surface output 2803 2871 2804 2872 n_out = 0 2805 2873 DO WHILE ( dosurf(1,n_out+1)(1:1) /= ' ' ) 2806 2874 2807 2875 n_out = n_out + 1 2808 2876 trimvar = TRIM( dosurf(1,n_out) ) … … 2826 2894 surf_def_v(3)%us, & 2827 2895 surf_lsm_v(3)%us, & 2828 surf_usm_v(3)%us, n_out ) 2896 surf_usm_v(3)%us, n_out ) 2829 2897 2830 2898 CASE ( 'ts' ) … … 2844 2912 surf_def_v(3)%ts, & 2845 2913 surf_lsm_v(3)%ts, & 2846 surf_usm_v(3)%ts, n_out ) 2914 surf_usm_v(3)%ts, n_out ) 2847 2915 2848 2916 CASE ( 'qs' ) … … 2862 2930 surf_def_v(3)%qs, & 2863 2931 surf_lsm_v(3)%qs, & 2864 surf_usm_v(3)%qs, n_out ) 2932 surf_usm_v(3)%qs, n_out ) 2865 2933 2866 2934 CASE ( 'ss' ) … … 2880 2948 surf_def_v(3)%ss, & 2881 2949 surf_lsm_v(3)%ss, & 2882 surf_usm_v(3)%ss, n_out ) 2950 surf_usm_v(3)%ss, n_out ) 2883 2951 2884 2952 CASE ( 'qcs' ) … … 2898 2966 surf_def_v(3)%qcs, & 2899 2967 surf_lsm_v(3)%qcs, & 2900 surf_usm_v(3)%qcs, n_out ) 2968 surf_usm_v(3)%qcs, n_out ) 2901 2969 2902 2970 CASE ( 'ncs' ) … … 2916 2984 surf_def_v(3)%ncs, & 2917 2985 surf_lsm_v(3)%ncs, & 2918 surf_usm_v(3)%ncs, n_out ) 2986 surf_usm_v(3)%ncs, n_out ) 2919 2987 2920 2988 CASE ( 'qrs' ) … … 2934 3002 surf_def_v(3)%qrs, & 2935 3003 surf_lsm_v(3)%qrs, & 2936 surf_usm_v(3)%qrs, n_out ) 3004 surf_usm_v(3)%qrs, n_out ) 2937 3005 2938 3006 CASE ( 'nrs' ) … … 2952 3020 surf_def_v(3)%nrs, & 2953 3021 surf_lsm_v(3)%nrs, & 2954 surf_usm_v(3)%nrs, n_out ) 3022 surf_usm_v(3)%nrs, n_out ) 2955 3023 2956 3024 CASE ( 'ol' ) … … 2970 3038 surf_def_v(3)%ol, & 2971 3039 surf_lsm_v(3)%ol, & 2972 surf_usm_v(3)%ol, n_out ) 3040 surf_usm_v(3)%ol, n_out ) 2973 3041 2974 3042 CASE ( 'z0' ) … … 2988 3056 surf_def_v(3)%z0, & 2989 3057 surf_lsm_v(3)%z0, & 2990 surf_usm_v(3)%z0, n_out ) 3058 surf_usm_v(3)%z0, n_out ) 2991 3059 2992 3060 CASE ( 'z0h' ) … … 3006 3074 surf_def_v(3)%z0h, & 3007 3075 surf_lsm_v(3)%z0h, & 3008 surf_usm_v(3)%z0h, n_out ) 3076 surf_usm_v(3)%z0h, n_out ) 3009 3077 3010 3078 CASE ( 'z0q' ) … … 3024 3092 surf_def_v(3)%z0q, & 3025 3093 surf_lsm_v(3)%z0q, & 3026 surf_usm_v(3)%z0q, n_out ) 3027 3028 CASE ( 'theta1' ) 3094 surf_usm_v(3)%z0q, n_out ) 3095 3096 CASE ( 'theta1' ) 3029 3097 CALL surface_data_output_sum_up( surf_def_h(0)%pt1, & 3030 3098 surf_def_h(1)%pt1, & … … 3042 3110 surf_def_v(3)%pt1, & 3043 3111 surf_lsm_v(3)%pt1, & 3044 surf_usm_v(3)%pt1, n_out ) 3045 3046 CASE ( 'qv1' ) 3112 surf_usm_v(3)%pt1, n_out ) 3113 3114 CASE ( 'qv1' ) 3047 3115 CALL surface_data_output_sum_up( surf_def_h(0)%qv1, & 3048 3116 surf_def_h(1)%qv1, & … … 3060 3128 surf_def_v(3)%qv1, & 3061 3129 surf_lsm_v(3)%qv1, & 3062 surf_usm_v(3)%qv1, n_out ) 3130 surf_usm_v(3)%qv1, n_out ) 3063 3131 3064 3132 CASE ( 'thetav1' ) … … 3078 3146 surf_def_v(3)%vpt1, & 3079 3147 surf_lsm_v(3)%vpt1, & 3080 surf_usm_v(3)%vpt1, n_out ) 3148 surf_usm_v(3)%vpt1, n_out ) 3081 3149 3082 3150 CASE ( 'usws' ) … … 3096 3164 surf_def_v(3)%usws, & 3097 3165 surf_lsm_v(3)%usws, & 3098 surf_usm_v(3)%usws, n_out ) 3166 surf_usm_v(3)%usws, n_out ) 3099 3167 3100 3168 CASE ( 'vsws' ) … … 3114 3182 surf_def_v(3)%vsws, & 3115 3183 surf_lsm_v(3)%vsws, & 3116 surf_usm_v(3)%vsws, n_out ) 3184 surf_usm_v(3)%vsws, n_out ) 3117 3185 3118 3186 CASE ( 'shf' ) … … 3120 3188 surf_def_h(1)%shf, & 3121 3189 surf_lsm_h%shf, & 3122 surf_usm_h%shf, & 3190 surf_usm_h%shf, & 3123 3191 surf_def_v(0)%shf, & 3124 3192 surf_lsm_v(0)%shf, & … … 3150 3218 surf_def_v(3)%qsws, & 3151 3219 surf_lsm_v(3)%qsws, & 3152 surf_usm_v(3)%qsws, n_out ) 3220 surf_usm_v(3)%qsws, n_out ) 3153 3221 3154 3222 CASE ( 'ssws' ) … … 3168 3236 surf_def_v(3)%ssws, & 3169 3237 surf_lsm_v(3)%ssws, & 3170 surf_usm_v(3)%ssws, n_out ) 3238 surf_usm_v(3)%ssws, n_out ) 3171 3239 3172 3240 CASE ( 'qcsws' ) … … 3186 3254 surf_def_v(3)%qcsws, & 3187 3255 surf_lsm_v(3)%qcsws, & 3188 surf_usm_v(3)%qcsws, n_out ) 3256 surf_usm_v(3)%qcsws, n_out ) 3189 3257 3190 3258 CASE ( 'ncsws' ) … … 3204 3272 surf_def_v(3)%ncsws, & 3205 3273 surf_lsm_v(3)%ncsws, & 3206 surf_usm_v(3)%ncsws, n_out ) 3274 surf_usm_v(3)%ncsws, n_out ) 3207 3275 3208 3276 CASE ( 'qrsws' ) … … 3222 3290 surf_def_v(3)%qrsws, & 3223 3291 surf_lsm_v(3)%qrsws, & 3224 surf_usm_v(3)%qrsws, n_out ) 3292 surf_usm_v(3)%qrsws, n_out ) 3225 3293 3226 3294 CASE ( 'nrsws' ) … … 3240 3308 surf_def_v(3)%nrsws, & 3241 3309 surf_lsm_v(3)%nrsws, & 3242 surf_usm_v(3)%nrsws, n_out ) 3310 surf_usm_v(3)%nrsws, n_out ) 3243 3311 3244 3312 CASE ( 'sasws' ) … … 3258 3326 surf_def_v(3)%sasws, & 3259 3327 surf_lsm_v(3)%sasws, & 3260 surf_usm_v(3)%sasws, n_out ) 3328 surf_usm_v(3)%sasws, n_out ) 3261 3329 3262 3330 CASE ( 'q_surface' ) … … 3276 3344 surf_def_v(3)%q_surface, & 3277 3345 surf_lsm_v(3)%q_surface, & 3278 surf_usm_v(3)%q_surface, n_out ) 3279 3346 surf_usm_v(3)%q_surface, n_out ) 3347 3280 3348 3281 3349 CASE ( 'theta_surface' ) … … 3295 3363 surf_def_v(3)%pt_surface, & 3296 3364 surf_lsm_v(3)%pt_surface, & 3297 surf_usm_v(3)%pt_surface, n_out ) 3365 surf_usm_v(3)%pt_surface, n_out ) 3298 3366 3299 3367 CASE ( 'thetav_surface' ) … … 3313 3381 surf_def_v(3)%vpt_surface, & 3314 3382 surf_lsm_v(3)%vpt_surface, & 3315 surf_usm_v(3)%vpt_surface, n_out ) 3383 surf_usm_v(3)%vpt_surface, n_out ) 3316 3384 3317 3385 CASE ( 'rad_net' ) … … 3331 3399 surf_def_v(3)%rad_net, & 3332 3400 surf_lsm_v(3)%rad_net, & 3333 surf_usm_v(3)%rad_net, n_out ) 3401 surf_usm_v(3)%rad_net, n_out ) 3334 3402 3335 3403 CASE ( 'rad_lw_in' ) … … 3349 3417 surf_def_v(3)%rad_lw_in, & 3350 3418 surf_lsm_v(3)%rad_lw_in, & 3351 surf_usm_v(3)%rad_lw_in, n_out ) 3419 surf_usm_v(3)%rad_lw_in, n_out ) 3352 3420 3353 3421 CASE ( 'rad_lw_out' ) … … 3367 3435 surf_def_v(3)%rad_lw_out, & 3368 3436 surf_lsm_v(3)%rad_lw_out, & 3369 surf_usm_v(3)%rad_lw_out, n_out ) 3437 surf_usm_v(3)%rad_lw_out, n_out ) 3370 3438 3371 3439 CASE ( 'rad_sw_in' ) … … 3385 3453 surf_def_v(3)%rad_sw_in, & 3386 3454 surf_lsm_v(3)%rad_sw_in, & 3387 surf_usm_v(3)%rad_sw_in, n_out ) 3455 surf_usm_v(3)%rad_sw_in, n_out ) 3388 3456 3389 3457 CASE ( 'rad_sw_out' ) … … 3403 3471 surf_def_v(3)%rad_sw_out, & 3404 3472 surf_lsm_v(3)%rad_sw_out, & 3405 surf_usm_v(3)%rad_sw_out, n_out ) 3473 surf_usm_v(3)%rad_sw_out, n_out ) 3406 3474 3407 3475 CASE ( 'ghf' ) … … 3421 3489 surf_def_v(3)%ghf, & 3422 3490 surf_lsm_v(3)%ghf, & 3423 surf_usm_v(3)%ghf, n_out ) 3424 3425 CASE ( 'r_a' ) 3491 surf_usm_v(3)%ghf, n_out ) 3492 3493 CASE ( 'r_a' ) 3426 3494 CALL surface_data_output_sum_up( surf_def_h(0)%r_a, & 3427 3495 surf_def_h(1)%r_a, & … … 3439 3507 surf_def_v(3)%r_a, & 3440 3508 surf_lsm_v(3)%r_a, & 3441 surf_usm_v(3)%r_a, n_out ) 3442 3443 CASE ( 'r_soil' ) 3509 surf_usm_v(3)%r_a, n_out ) 3510 3511 CASE ( 'r_soil' ) 3444 3512 CALL surface_data_output_sum_up( surf_def_h(0)%r_soil, & 3445 3513 surf_def_h(1)%r_soil, & … … 3457 3525 surf_def_v(3)%r_soil, & 3458 3526 surf_lsm_v(3)%r_soil, & 3459 surf_usm_v(3)%r_soil, n_out ) 3527 surf_usm_v(3)%r_soil, n_out ) 3460 3528 3461 3529 CASE ( 'r_canopy' ) … … 3475 3543 surf_def_v(3)%r_canopy, & 3476 3544 surf_lsm_v(3)%r_canopy, & 3477 surf_usm_v(3)%r_canopy, n_out ) 3545 surf_usm_v(3)%r_canopy, n_out ) 3478 3546 3479 3547 CASE ( 'r_s' ) … … 3493 3561 surf_def_v(3)%r_s, & 3494 3562 surf_lsm_v(3)%r_s, & 3495 surf_usm_v(3)%r_s, n_out ) 3563 surf_usm_v(3)%r_s, n_out ) 3496 3564 3497 3565 … … 3548 3616 surf_lsm_v(3)%rad_sw_ref, & 3549 3617 surf_usm_v(3)%rad_sw_ref, n_out ) 3550 3618 3551 3619 CASE ( 'rad_sw_res' ) 3552 3620 CALL surface_data_output_sum_up( surf_def_h(0)%rad_sw_res, & … … 3620 3688 surf_lsm_v(3)%rad_lw_res, & 3621 3689 surf_usm_v(3)%rad_lw_res, n_out ) 3622 3690 3623 3691 CASE ( 'uvw1' ) 3624 3692 CALL surface_data_output_sum_up( surf_def_h(0)%uvw_abs, & … … 3637 3705 surf_def_v(3)%uvw_abs, & 3638 3706 surf_lsm_v(3)%uvw_abs, & 3639 surf_usm_v(3)%uvw_abs, n_out ) 3707 surf_usm_v(3)%uvw_abs, n_out ) 3640 3708 3641 3709 END SELECT 3642 3710 ENDDO 3643 3644 3711 3712 3645 3713 END SUBROUTINE surface_data_output_averaging 3646 3714 3647 3715 !! 3648 3716 ! Description: 3649 3717 !  3650 !> Sumup the surface data for average output variables. 3651 !! 3718 !> Sumup the surface data for average output variables. 3719 !! 3652 3720 SUBROUTINE surface_data_output_sum_up( var_def_h0, var_def_h1, & 3653 3721 var_lsm_h, var_usm_h, & … … 3656 3724 var_def_v2, var_lsm_v2, var_usm_v2, & 3657 3725 var_def_v3, var_lsm_v3, var_usm_v3, n_out ) 3658 3726 3659 3727 IMPLICIT NONE 3660 3728 … … 3662 3730 INTEGER(iwp) :: n_out !< index for output variable 3663 3731 INTEGER(iwp) :: n_surf !< running index for surface elements 3664 3732 3665 3733 REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(IN) :: var_def_h0 !< output variable at upwardfacing defaulttype surfaces 3666 3734 REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(IN) :: var_def_h1 !< output variable at downwardfacing defaulttype surfaces … … 3679 3747 REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(IN) :: var_usm_v2 !< output variable at eastwardfacing urbantype surfaces 3680 3748 REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(IN) :: var_usm_v3 !< output variable at westwardfacing urbantype surfaces 3681 3682 ! 3683 ! Set counter variable to zero before the variable is written to 3684 ! the output array. 3749 3750 ! 3751 ! Set counter variable to zero before the variable is written to 3752 ! the output array. 3685 3753 n_surf = 0 3686 3754 3687 3755 ! 3688 3756 ! Write the horizontal surfaces. 3689 ! Before each the variable is written to the output data structure, first 3690 ! check if the variable for the respective surface type is defined. 3757 ! Before each the variable is written to the output data structure, first 3758 ! check if the variable for the respective surface type is defined. 3691 3759 ! If a variable is not defined, skip the block and increment the counter 3692 ! variable by the number of surface elements of this type. Usually this 3760 ! variable by the number of surface elements of this type. Usually this 3693 3761 ! is zere, however, there might be the situation that e.g. urban surfaces 3694 3762 ! are defined but the respective variable is not allocated for this surface 3695 ! type. To write the data on the exact position, increment the counter. 3763 ! type. To write the data on the exact position, increment the counter. 3696 3764 IF ( ALLOCATED( var_def_h0 ) ) THEN 3697 DO m = 1, surf_def_h(0)%ns 3765 DO m = 1, surf_def_h(0)%ns 3698 3766 n_surf = n_surf + 1 3699 3767 surfaces%var_av(n_surf,n_out) = surfaces%var_av(n_surf,n_out) & … … 3701 3769 ENDDO 3702 3770 ELSE 3703 n_surf = n_surf + surf_def_h(0)%ns 3771 n_surf = n_surf + surf_def_h(0)%ns 3704 3772 ENDIF 3705 3773 IF ( ALLOCATED( var_def_h1 ) ) THEN 3706 DO m = 1, surf_def_h(1)%ns 3774 DO m = 1, surf_def_h(1)%ns 3707 3775 n_surf = n_surf + 1 3708 3776 surfaces%var_av(n_surf,n_out) = surfaces%var_av(n_surf,n_out) & … … 3713 3781 ENDIF 3714 3782 IF ( ALLOCATED( var_lsm_h ) ) THEN 3715 DO m = 1, surf_lsm_h%ns 3783 DO m = 1, surf_lsm_h%ns 3716 3784 n_surf = n_surf + 1 3717 3785 surfaces%var_av(n_surf,n_out) = surfaces%var_av(n_surf,n_out) & … … 3722 3790 ENDIF 3723 3791 IF ( ALLOCATED( var_usm_h ) ) THEN 3724 DO m = 1, surf_usm_h%ns 3792 DO m = 1, surf_usm_h%ns 3725 3793 n_surf = n_surf + 1 3726 3794 surfaces%var_av(n_surf,n_out) = surfaces%var_av(n_surf,n_out) & … … 3733 3801 ! Write northwardfacing 3734 3802 IF ( ALLOCATED( var_def_v0 ) ) THEN 3735 DO m = 1, surf_def_v(0)%ns 3803 DO m = 1, surf_def_v(0)%ns 3736 3804 n_surf = n_surf + 1 3737 3805 surfaces%var_av(n_surf,n_out) = surfaces%var_av(n_surf,n_out) & … … 3742 3810 ENDIF 3743 3811 IF ( ALLOCATED( var_lsm_v0 ) ) THEN 3744 DO m = 1, surf_lsm_v(0)%ns 3812 DO m = 1, surf_lsm_v(0)%ns 3745 3813 n_surf = n_surf + 1 3746 3814 surfaces%var_av(n_surf,n_out) = surfaces%var_av(n_surf,n_out) & … … 3751 3819 ENDIF 3752 3820 IF ( ALLOCATED( var_usm_v0 ) ) THEN 3753 DO m = 1, surf_usm_v(0)%ns 3821 DO m = 1, surf_usm_v(0)%ns 3754 3822 n_surf = n_surf + 1 3755 3823 surfaces%var_av(n_surf,n_out) = surfaces%var_av(n_surf,n_out) & … … 3762 3830 ! Write southwardfacing 3763 3831 IF ( ALLOCATED( var_def_v1 ) ) THEN 3764 DO m = 1, surf_def_v(1)%ns 3832 DO m = 1, surf_def_v(1)%ns 3765 3833 n_surf = n_surf + 1 3766 3834 surfaces%var_av(n_surf,n_out) = surfaces%var_av(n_surf,n_out) & … … 3771 3839 ENDIF 3772 3840 IF ( ALLOCATED( var_lsm_v1 ) ) THEN 3773 DO m = 1, surf_lsm_v(1)%ns 3841 DO m = 1, surf_lsm_v(1)%ns 3774 3842 n_surf = n_surf + 1 3775 3843 surfaces%var_av(n_surf,n_out) = surfaces%var_av(n_surf,n_out) & … … 3780 3848 ENDIF 3781 3849 IF ( ALLOCATED( var_usm_v1 ) ) THEN 3782 DO m = 1, surf_usm_v(1)%ns 3850 DO m = 1, surf_usm_v(1)%ns 3783 3851 n_surf = n_surf + 1 3784 3852 surfaces%var_av(n_surf,n_out) = surfaces%var_av(n_surf,n_out) & … … 3791 3859 ! Write eastwardfacing 3792 3860 IF ( ALLOCATED( var_def_v2 ) ) THEN 3793 DO m = 1, surf_def_v(2)%ns 3861 DO m = 1, surf_def_v(2)%ns 3794 3862 n_surf = n_surf + 1 3795 3863 surfaces%var_av(n_surf,n_out) = surfaces%var_av(n_surf,n_out) & … … 3800 3868 ENDIF 3801 3869 IF ( ALLOCATED( var_lsm_v2 ) ) THEN 3802 DO m = 1, surf_lsm_v(2)%ns 3870 DO m = 1, surf_lsm_v(2)%ns 3803 3871 n_surf = n_surf + 1 3804 3872 surfaces%var_av(n_surf,n_out) = surfaces%var_av(n_surf,n_out) & … … 3809 3877 ENDIF 3810 3878 IF ( ALLOCATED( var_usm_v2 ) ) THEN 3811 DO m = 1, surf_usm_v(2)%ns 3879 DO m = 1, surf_usm_v(2)%ns 3812 3880 n_surf = n_surf + 1 3813 3881 surfaces%var_av(n_surf,n_out) = surfaces%var_av(n_surf,n_out) & … … 3820 3888 ! Write westwardfacing 3821 3889 IF ( ALLOCATED( var_def_v3 ) ) THEN 3822 DO m = 1, surf_def_v(3)%ns 3890 DO m = 1, surf_def_v(3)%ns 3823 3891 n_surf = n_surf + 1 3824 3892 surfaces%var_av(n_surf,n_out) = surfaces%var_av(n_surf,n_out) & … … 3829 3897 ENDIF 3830 3898 IF ( ALLOCATED( var_lsm_v3 ) ) THEN 3831 DO m = 1, surf_lsm_v(3)%ns 3899 DO m = 1, surf_lsm_v(3)%ns 3832 3900 n_surf = n_surf + 1 3833 3901 surfaces%var_av(n_surf,n_out) = surfaces%var_av(n_surf,n_out) & … … 3838 3906 ENDIF 3839 3907 IF ( ALLOCATED( var_usm_v3 ) ) THEN 3840 DO m = 1, surf_usm_v(3)%ns 3908 DO m = 1, surf_usm_v(3)%ns 3841 3909 n_surf = n_surf + 1 3842 3910 surfaces%var_av(n_surf,n_out) = surfaces%var_av(n_surf,n_out) & … … 3846 3914 n_surf = n_surf + surf_usm_v(3)%ns 3847 3915 ENDIF 3848 3916 3849 3917 END SUBROUTINE surface_data_output_sum_up 3850 3918 3851 3919 !! 3852 3920 ! Description: … … 3859 3927 var_def_v1, var_lsm_v1, var_usm_v1, & 3860 3928 var_def_v2, var_lsm_v2, var_usm_v2, & 3861 var_def_v3, var_lsm_v3, var_usm_v3 ) 3862 3929 var_def_v3, var_lsm_v3, var_usm_v3 ) 3930 3863 3931 IMPLICIT NONE 3864 3932 3865 3933 INTEGER(iwp) :: m !< running index for surface elements 3866 3934 INTEGER(iwp) :: n_surf !< running index for surface elements 3867 3935 3868 3936 REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(IN) :: var_def_h0 !< output variable at upwardfacing defaulttype surfaces 3869 3937 REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(IN) :: var_def_h1 !< output variable at downwardfacing defaulttype surfaces … … 3882 3950 REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(IN) :: var_usm_v2 !< output variable at eastwardfacing urbantype surfaces 3883 3951 REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(IN) :: var_usm_v3 !< output variable at westwardfacing urbantype surfaces 3884 3885 ! 3886 ! Set counter variable to zero before the variable is written to 3887 ! the output array. 3952 3953 ! 3954 ! Set counter variable to zero before the variable is written to 3955 ! the output array. 3888 3956 n_surf = 0 3889 3957 3890 3958 ! 3891 3959 ! Write the horizontal surfaces. 3892 ! Before each the variable is written to the output data structure, first 3893 ! check if the variable for the respective surface type is defined. 3960 ! Before each the variable is written to the output data structure, first 3961 ! check if the variable for the respective surface type is defined. 3894 3962 ! If a variable is not defined, skip the block and increment the counter 3895 ! variable by the number of surface elements of this type. Usually this 3963 ! variable by the number of surface elements of this type. Usually this 3896 3964 ! is zere, however, there might be the situation that e.g. urban surfaces 3897 3965 ! are defined but the respective variable is not allocated for this surface 3898 ! type. To write the data on the exact position, increment the counter. 3966 ! type. To write the data on the exact position, increment the counter. 3899 3967 IF ( ALLOCATED( var_def_h0 ) ) THEN 3900 DO m = 1, surf_def_h(0)%ns 3968 DO m = 1, surf_def_h(0)%ns 3901 3969 n_surf = n_surf + 1 3902 3970 surfaces%var_out(n_surf) = var_def_h0(m) 3903 3971 ENDDO 3904 3972 ELSE 3905 n_surf = n_surf + surf_def_h(0)%ns 3973 n_surf = n_surf + surf_def_h(0)%ns 3906 3974 ENDIF 3907 3975 IF ( ALLOCATED( var_def_h1 ) ) THEN 3908 DO m = 1, surf_def_h(1)%ns 3976 DO m = 1, surf_def_h(1)%ns 3909 3977 n_surf = n_surf + 1 3910 3978 surfaces%var_out(n_surf) = var_def_h1(m) … … 3914 3982 ENDIF 3915 3983 IF ( ALLOCATED( var_lsm_h ) ) THEN 3916 DO m = 1, surf_lsm_h%ns 3984 DO m = 1, surf_lsm_h%ns 3917 3985 n_surf = n_surf + 1 3918 3986 surfaces%var_out(n_surf) = var_lsm_h(m) … … 3922 3990 ENDIF 3923 3991 IF ( ALLOCATED( var_usm_h ) ) THEN 3924 DO m = 1, surf_usm_h%ns 3992 DO m = 1, surf_usm_h%ns 3925 3993 n_surf = n_surf + 1 3926 3994 surfaces%var_out(n_surf) = var_usm_h(m) … … 3932 4000 ! Write northwardfacing 3933 4001 IF ( ALLOCATED( var_def_v0 ) ) THEN 3934 DO m = 1, surf_def_v(0)%ns 4002 DO m = 1, surf_def_v(0)%ns 3935 4003 n_surf = n_surf + 1 3936 4004 surfaces%var_out(n_surf) = var_def_v0(m) … … 3940 4008 ENDIF 3941 4009 IF ( ALLOCATED( var_lsm_v0 ) ) THEN 3942 DO m = 1, surf_lsm_v(0)%ns 4010 DO m = 1, surf_lsm_v(0)%ns 3943 4011 n_surf = n_surf + 1 3944 4012 surfaces%var_out(n_surf) = var_lsm_v0(m) … … 3948 4016 ENDIF 3949 4017 IF ( ALLOCATED( var_usm_v0 ) ) THEN 3950 DO m = 1, surf_usm_v(0)%ns 4018 DO m = 1, surf_usm_v(0)%ns 3951 4019 n_surf = n_surf + 1 3952 4020 surfaces%var_out(n_surf) = var_usm_v0(m) … … 3958 4026 ! Write southwardfacing 3959 4027 IF ( ALLOCATED( var_def_v1 ) ) THEN 3960 DO m = 1, surf_def_v(1)%ns 4028 DO m = 1, surf_def_v(1)%ns 3961 4029 n_surf = n_surf + 1 3962 4030 surfaces%var_out(n_surf) = var_def_v1(m) … … 3966 4034 ENDIF 3967 4035 IF ( ALLOCATED( var_lsm_v1 ) ) THEN 3968 DO m = 1, surf_lsm_v(1)%ns 4036 DO m = 1, surf_lsm_v(1)%ns 3969 4037 n_surf = n_surf + 1 3970 4038 surfaces%var_out(n_surf) = var_lsm_v1(m) … … 3974 4042 ENDIF 3975 4043 IF ( ALLOCATED( var_usm_v1 ) ) THEN 3976 DO m = 1, surf_usm_v(1)%ns 4044 DO m = 1, surf_usm_v(1)%ns 3977 4045 n_surf = n_surf + 1 3978 4046 surfaces%var_out(n_surf) = var_usm_v1(m) … … 3984 4052 ! Write eastwardfacing 3985 4053 IF ( ALLOCATED( var_def_v2 ) ) THEN 3986 DO m = 1, surf_def_v(2)%ns 4054 DO m = 1, surf_def_v(2)%ns 3987 4055 n_surf = n_surf + 1 3988 4056 surfaces%var_out(n_surf) = var_def_v2(m) … … 3992 4060 ENDIF 3993 4061 IF ( ALLOCATED( var_lsm_v2 ) ) THEN 3994 DO m = 1, surf_lsm_v(2)%ns 4062 DO m = 1, surf_lsm_v(2)%ns 3995 4063 n_surf = n_surf + 1 3996 4064 surfaces%var_out(n_surf) = var_lsm_v2(m) … … 4000 4068 ENDIF 4001 4069 IF ( ALLOCATED( var_usm_v2 ) ) THEN 4002 DO m = 1, surf_usm_v(2)%ns 4070 DO m = 1, surf_usm_v(2)%ns 4003 4071 n_surf = n_surf + 1 4004 4072 surfaces%var_out(n_surf) = var_usm_v2(m) … … 4010 4078 ! Write westwardfacing 4011 4079 IF ( ALLOCATED( var_def_v3 ) ) THEN 4012 DO m = 1, surf_def_v(3)%ns 4080 DO m = 1, surf_def_v(3)%ns 4013 4081 n_surf = n_surf + 1 4014 4082 surfaces%var_out(n_surf) = var_def_v3(m) … … 4018 4086 ENDIF 4019 4087 IF ( ALLOCATED( var_lsm_v3 ) ) THEN 4020 DO m = 1, surf_lsm_v(3)%ns 4088 DO m = 1, surf_lsm_v(3)%ns 4021 4089 n_surf = n_surf + 1 4022 4090 surfaces%var_out(n_surf) = var_lsm_v3(m) … … 4026 4094 ENDIF 4027 4095 IF ( ALLOCATED( var_usm_v3 ) ) THEN 4028 DO m = 1, surf_usm_v(3)%ns 4096 DO m = 1, surf_usm_v(3)%ns 4029 4097 n_surf = n_surf + 1 4030 4098 surfaces%var_out(n_surf) = var_usm_v3(m) … … 4033 4101 n_surf = n_surf + surf_usm_v(3)%ns 4034 4102 ENDIF 4035 4103 4036 4104 END SUBROUTINE surface_data_output_collect 4037 4105 4038 4106 !! 4039 4107 ! Description: … … 4045 4113 IMPLICIT NONE 4046 4114 4047 CHARACTER (LEN=80) :: line !< dummy string that contains the current line of the parameter file 4048 4115 CHARACTER (LEN=80) :: line !< dummy string that contains the current line of the parameter file 4116 4049 4117 4050 4118 NAMELIST /surface_data_output_parameters/ & … … 4053 4121 skip_time_dosurf, skip_time_dosurf_av, & 4054 4122 to_netcdf, to_vtk 4055 4123 4056 4124 line = ' ' 4057 4125 4058 4126 ! 4059 4127 ! Try to find the namelist … … 4075 4143 10 BACKSPACE( 11 ) 4076 4144 READ( 11 , '(A)') line 4077 CALL parin_fail_message( 'surface_data_output_parameters', line ) 4078 4145 CALL parin_fail_message( 'surface_data_output_parameters', line ) 4146 4079 4147 14 CONTINUE 4080 4148 4081 4149 4082 4150 END SUBROUTINE surface_data_output_parin 4083 4084 4151 4152 4085 4153 !! 4086 4154 ! Description: 4087 4155 !  4088 !> Check the input parameters for consistency. Further preprocess the given 4156 !> Check the input parameters for consistency. Further preprocess the given 4089 4157 !> output variables, i.e. separate them into average and nonaverage output 4090 !> variables and map them onto internal output array. 4158 !> variables and map them onto internal output array. 4091 4159 !! 4092 4160 SUBROUTINE surface_data_output_check_parameters … … 4095 4163 ONLY: averaging_interval, dt_data_output, initializing_actions, & 4096 4164 message_string 4097 4165 4098 4166 USE pegrid, & 4099 4167 ONLY: numprocs_previous_run … … 4125 4193 'PA0536', 1, 2, 0, 6, 0 ) 4126 4194 ENDIF 4127 4195 4128 4196 #if ! defined( __netcdf4_parallel ) 4129 4197 ! … … 4147 4215 ENDIF 4148 4216 ! 4149 ! Count number of output variables and separate output strings for 4217 ! Count number of output variables and separate output strings for 4150 4218 ! average and nonaverage output variables. 4151 4219 n_out = 0 … … 4194 4262 CASE ( 'usws', 'vsws' ) 4195 4263 unit = 'm2/s2' 4196 4264 4197 4265 CASE ( 'qcsws', 'ncsws', 'qrsws', 'nrsws', 'sasws' ) 4198 4266 … … 4238 4306 4239 4307 END SUBROUTINE surface_data_output_check_parameters 4240 4241 4308 4309 4242 4310 !! 4243 4311 ! Description: … … 4246 4314 !! 4247 4315 SUBROUTINE surface_data_output_last_action( av ) 4248 4316 4249 4317 USE control_parameters, & 4250 4318 ONLY: io_blocks, io_group … … 4254 4322 4255 4323 IMPLICIT NONE 4256 4324 4257 4325 INTEGER(iwp) :: av !< id indicating average or nonaverage data output 4258 4326 INTEGER(iwp) :: i !< loop index … … 4276 4344 CALL MPI_BARRIER( comm2d, ierr ) 4277 4345 #endif 4278 ENDDO 4346 ENDDO 4279 4347 ENDIF 4280 4348 4281 END SUBROUTINE surface_data_output_last_action 4282 4349 END SUBROUTINE surface_data_output_last_action 4350 4283 4351 !! 4284 4352 ! Description: … … 4291 4359 USE control_parameters, & 4292 4360 ONLY: length, restart_string 4293 4361 4294 4362 IMPLICIT NONE 4295 4363 4296 4364 LOGICAL, INTENT(OUT) :: found !< flag indicating if variable was found 4297 4365 4298 4366 found = .TRUE. 4299 4367 4300 4368 SELECT CASE ( restart_string(1:length) ) 4301 4369 4302 4370 CASE ( 'average_count_surf' ) 4303 4371 READ ( 13 ) average_count_surf 4304 4372 4305 4373 CASE DEFAULT 4306 4374 4307 4375 found = .FALSE. 4308 4376 4309 4377 END SELECT 4310 4378 4311 4379 4312 4380 END SUBROUTINE surface_data_output_rrd_global 4313 4381 4314 4382 !! 4315 4383 ! Description: … … 4324 4392 USE control_parameters, & 4325 4393 ONLY: length, restart_string 4326 4394 4327 4395 IMPLICIT NONE 4328 4396 … … 4332 4400 INTEGER(iwp) :: ns_h_on_file_usm !< number of horizontal surface elements (urban type) on file 4333 4401 INTEGER(iwp) :: nxlc !< index of left boundary on current subdomain 4334 INTEGER(iwp) :: nxlf !< index of left boundary on former subdomain 4335 INTEGER(iwp) :: nxl_on_file !< index of left boundary on former local domain 4402 INTEGER(iwp) :: nxlf !< index of left boundary on former subdomain 4403 INTEGER(iwp) :: nxl_on_file !< index of left boundary on former local domain 4336 4404 INTEGER(iwp) :: nxrc !< index of right boundary on current subdomain 4337 4405 INTEGER(iwp) :: nxrf !< index of right boundary on former subdomain 4338 INTEGER(iwp) :: nxr_on_file !< index of right boundary on former local domain 4406 INTEGER(iwp) :: nxr_on_file !< index of right boundary on former local domain 4339 4407 INTEGER(iwp) :: nync !< index of north boundary on current subdomain 4340 4408 INTEGER(iwp) :: nynf !< index of north boundary on former subdomain 4341 INTEGER(iwp) :: nyn_on_file !< index of north boundary on former local domain 4342 INTEGER(iwp) :: nysc !< index of south boundary on current subdomain 4409 INTEGER(iwp) :: nyn_on_file !< index of north boundary on former local domain 4410 INTEGER(iwp) :: nysc !< index of south boundary on current subdomain 4343 4411 INTEGER(iwp) :: nysf !< index of south boundary on former subdomain 4344 4412 INTEGER(iwp) :: nys_on_file !< index of south boundary on former local domain … … 4354 4422 CASE ( 'surfaces%var_av' ) 4355 4423 IF ( k == 1 ) READ ( 13 ) surfaces%var_av 4356 4424 4357 4425 CASE DEFAULT 4358 4426 … … 4363 4431 4364 4432 END SUBROUTINE surface_data_output_rrd_local 4365 4433 4366 4434 !! 4367 4435 ! Description: … … 4377 4445 4378 4446 END SUBROUTINE surface_data_output_wrd_global 4379 4447 4380 4448 !! 4381 4449 ! Description:
Note: See TracChangeset
for help on using the changeset viewer.