Changeset 3494 for palm/trunk/SOURCE/surface_output_mod.f90
- Timestamp:
- Nov 6, 2018 2:51:27 PM (6 years ago)
- Location:
- palm/trunk/SOURCE
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
surface_output_mod.f90 (modified) (27 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE
- Property svn:mergeinfo changed
/palm/branches/mosaik_M2 (added) merged: 2360,3363,3434,3437,3471
- Property svn:mergeinfo changed
-
palm/trunk/SOURCE/surface_output_mod.f90
r3483 r3494 25 25 ! ----------------- 26 26 ! $Id$ 27 ! Bugfix in gathering surface data from different types and orientation. 28 ! Output of total number of surfaces and vertices added. 29 ! String length is output, for more convinient post-processing. 30 ! Last actions added. 31 ! 32 ! 3483 2018-11-02 14:19:26Z raasch 27 33 ! bugfix: missed directives for MPI added 28 34 ! … … 30 36 ! Add output variables 31 37 ! Write data into binary file 38 ! 39 ! 3420 2018-10-24 17:30:08Z gronemeier 32 40 ! Initial implementation from Klaus Ketelsen and Matthias Suehring 33 41 ! 34 ! 3420 2018-10-24 17:30:08Z gronemeier 42 ! 43 ! Authors: 44 ! -------- 45 ! @author Klaus Ketelsen, Matthias Suehring 35 46 ! 36 47 ! Description: 37 48 ! ------------ 38 49 !> Generate output for surface data. 39 ! 50 !> 51 !> @todo Create namelist file for post-processing tool. 40 52 !------------------------------------------------------------------------------! 41 53 … … 46 58 USE arrays_3d, & 47 59 ONLY: zu, zw 60 48 61 USE control_parameters, & 49 62 ONLY: surface_data_output 63 50 64 USE grid_variables, & 51 65 ONLY: dx,dy 66 52 67 USE indices, & 53 68 ONLY: nxl, nxr, nys, nyn, nzb, nzt 69 70 USE pegrid 71 54 72 USE surface_mod, & 55 73 ONLY: surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, & … … 58 76 IMPLICIT NONE 59 77 60 TYPE surf_out !< data structure which contains all surfaces elements of all types on subdomain78 TYPE surf_out !< data structure which contains all surfaces elements of all types on subdomain 61 79 62 INTEGER(iwp) :: ns !< number of surface elements 63 INTEGER(iwp) :: npoints !< number of points / vertices which define a surface element 80 INTEGER(iwp) :: ns !< number of surface elements on subdomain 81 INTEGER(iwp) :: ns_total !< total number of surface elements 82 INTEGER(iwp) :: npoints !< number of points / vertices which define a surface element (on subdomain) 83 INTEGER(iwp) :: npoints_total !< total number of points / vertices which define a surface element 64 84 65 85 REAL(wp) :: fillvalue = -9999.9_wp !< fillvalue for surface elements which are not defined … … 111 131 END INTERFACE surface_output_init 112 132 133 INTERFACE surface_output_last_action 134 MODULE PROCEDURE surface_output_last_action 135 END INTERFACE surface_output_last_action 136 113 137 INTERFACE surface_output_parin 114 138 MODULE PROCEDURE surface_output_parin … … 119 143 PUBLIC surface_output, surface_output_averaging, & 120 144 surface_output_check_parameters, surface_output_init, & 121 surface_output_ parin145 surface_output_last_action, surface_output_parin 122 146 ! 123 147 !--Public variables … … 134 158 !------------------------------------------------------------------------------! 135 159 SUBROUTINE surface_output_init 160 136 161 IMPLICIT NONE 137 162 138 INTEGER(iwp) :: i !< grid index in x-direction, also running variable for counting non-average data output 139 INTEGER(iwp) :: ilen !< string length 140 INTEGER(iwp) :: j !< grid index in y-direction, also running variable for counting average data output 141 INTEGER(iwp) :: k !< grid index in z-direction 142 INTEGER(iwp) :: l !< running index for surface-element orientation 143 INTEGER(iwp) :: m !< running index for surface elements 144 INTEGER(iwp) :: n_out !< running index for number of output variables 163 INTEGER(iwp) :: i !< grid index in x-direction, also running variable for counting non-average data output 164 INTEGER(iwp) :: ilen !< string length 165 INTEGER(iwp) :: j !< grid index in y-direction, also running variable for counting average data output 166 INTEGER(iwp) :: k !< grid index in z-direction 167 INTEGER(iwp) :: l !< running index for surface-element orientation 168 INTEGER(iwp) :: m !< running index for surface elements 169 INTEGER(iwp) :: n_out !< running index for number of output variables 170 INTEGER(iwp) :: npg !< counter variable for all surface elements ( or polygons ) 171 INTEGER(iwp) :: point_index_count !< local counter variable for point index 145 172 146 147 INTEGER(iwp) :: npg !< counter variable for all surface elements ( or polygons ) 148 173 INTEGER(iwp), DIMENSION(0:numprocs-1) :: num_points_on_pe !< array which contains the number of points on all mpi ranks 149 174 INTEGER(iwp), ALLOCATABLE, DIMENSION(:,:,:) :: point_index !< dummy array used to check where the reference points for surface polygons are located 150 175 … … 187 212 !-- Check if the vertices that define the surface element are already 188 213 !-- defined, if not, increment the counter. 189 IF ( point_index(k,j,i) < 0 ) THEN214 IF ( point_index(k,j,i) < 0 ) THEN 190 215 surfaces%npoints = surfaces%npoints + 1 191 point_index(k,j,i) = surfaces%npoints 216 point_index(k,j,i) = surfaces%npoints - 1 192 217 ENDIF 193 IF ( point_index(k,j +1,i) < 0) THEN218 IF ( point_index(k,j,i+1) < 0 ) THEN 194 219 surfaces%npoints = surfaces%npoints + 1 195 point_index(k,j +1,i) = surfaces%npoints220 point_index(k,j,i+1) = surfaces%npoints - 1 196 221 ENDIF 197 IF ( point_index(k,j,i+1) < 0) THEN 222 IF ( point_index(k,j+1,i+1) < 0 ) THEN 223 surfaces%npoints = surfaces%npoints + 1 224 point_index(k,j+1,i+1) = surfaces%npoints - 1 225 ENDIF 226 IF ( point_index(k,j+1,i) < 0 ) THEN 198 227 surfaces%npoints = surfaces%npoints + 1 199 point_index(k,j,i+1) = surfaces%npoints 200 ENDIF 201 IF ( point_index(k,j+1,i+1) < 0) THEN 202 surfaces%npoints = surfaces%npoints + 1 203 point_index(k,j+1,i+1) = surfaces%npoints 204 ENDIF 228 point_index(k,j+1,i) = surfaces%npoints - 1 229 ENDIF 205 230 ENDDO 206 231 ENDDO … … 210 235 k = surf_lsm_h%k(m) + surf_lsm_h%koff 211 236 212 IF ( point_index(k,j,i) < 0 ) THEN237 IF ( point_index(k,j,i) < 0 ) THEN 213 238 surfaces%npoints = surfaces%npoints + 1 214 point_index(k,j,i) = surfaces%npoints 239 point_index(k,j,i) = surfaces%npoints - 1 215 240 ENDIF 216 IF ( point_index(k,j +1,i) < 0) THEN241 IF ( point_index(k,j,i+1) < 0 ) THEN 217 242 surfaces%npoints = surfaces%npoints + 1 218 point_index(k,j +1,i) = surfaces%npoints243 point_index(k,j,i+1) = surfaces%npoints - 1 219 244 ENDIF 220 IF ( point_index(k,j,i+1) < 0) THEN 245 IF ( point_index(k,j+1,i+1) < 0 ) THEN 246 surfaces%npoints = surfaces%npoints + 1 247 point_index(k,j+1,i+1) = surfaces%npoints - 1 248 ENDIF 249 IF ( point_index(k,j+1,i) < 0 ) THEN 221 250 surfaces%npoints = surfaces%npoints + 1 222 point_index(k,j,i+1) = surfaces%npoints 223 ENDIF 224 IF ( point_index(k,j+1,i+1) < 0) THEN 225 surfaces%npoints = surfaces%npoints + 1 226 point_index(k,j+1,i+1) = surfaces%npoints 227 ENDIF 251 point_index(k,j+1,i) = surfaces%npoints - 1 252 ENDIF 228 253 ENDDO 229 254 DO m = 1, surf_usm_h%ns … … 232 257 k = surf_usm_h%k(m) + surf_usm_h%koff 233 258 234 IF ( point_index(k,j,i) < 0 ) THEN259 IF ( point_index(k,j,i) < 0 ) THEN 235 260 surfaces%npoints = surfaces%npoints + 1 236 point_index(k,j,i) = surfaces%npoints 261 point_index(k,j,i) = surfaces%npoints - 1 237 262 ENDIF 238 IF ( point_index(k,j +1,i) < 0) THEN263 IF ( point_index(k,j,i+1) < 0 ) THEN 239 264 surfaces%npoints = surfaces%npoints + 1 240 point_index(k,j +1,i) = surfaces%npoints265 point_index(k,j,i+1) = surfaces%npoints - 1 241 266 ENDIF 242 IF ( point_index(k,j,i+1) < 0) THEN 267 IF ( point_index(k,j+1,i+1) < 0 ) THEN 268 surfaces%npoints = surfaces%npoints + 1 269 point_index(k,j+1,i+1) = surfaces%npoints - 1 270 ENDIF 271 IF ( point_index(k,j+1,i) < 0 ) THEN 243 272 surfaces%npoints = surfaces%npoints + 1 244 point_index(k,j,i+1) = surfaces%npoints 245 ENDIF 246 IF ( point_index(k,j+1,i+1) < 0) THEN 247 surfaces%npoints = surfaces%npoints + 1 248 point_index(k,j+1,i+1) = surfaces%npoints 273 point_index(k,j+1,i) = surfaces%npoints - 1 249 274 ENDIF 250 275 ENDDO … … 258 283 !-- identical to the reference j-index outside the grid box. 259 284 !-- Equivalent for east-facing surfaces and i-index. 260 i = surf_def_v(l)%i(m) + MERGE( surf_def_v(l)%ioff, 0, l == 2)261 j = surf_def_v(l)%j(m) + MERGE( surf_def_v(l)%joff, 0, l == 0)285 i = surf_def_v(l)%i(m) + MERGE( surf_def_v(l)%ioff, 0, l == 3 ) 286 j = surf_def_v(l)%j(m) + MERGE( surf_def_v(l)%joff, 0, l == 1 ) 262 287 k = surf_def_v(l)%k(m) + surf_def_v(l)%koff 263 288 ! 264 289 !-- Lower left /front vertex 265 IF ( point_index(k,j,i) < 0 ) THEN290 IF ( point_index(k,j,i) < 0 ) THEN 266 291 surfaces%npoints = surfaces%npoints + 1 267 point_index(k,j,i) = surfaces%npoints 292 point_index(k,j,i) = surfaces%npoints - 1 268 293 ENDIF 269 294 ! 295 !-- Upper / lower right index for north- and south-facing surfaces 296 IF ( l == 0 .OR. l == 1 ) THEN 297 IF ( point_index(k,j,i+1) < 0 ) THEN 298 surfaces%npoints = surfaces%npoints + 1 299 point_index(k,j,i+1) = surfaces%npoints - 1 300 ENDIF 301 IF ( point_index(k+1,j,i+1) < 0 ) THEN 302 surfaces%npoints = surfaces%npoints + 1 303 point_index(k+1,j,i+1) = surfaces%npoints - 1 304 ENDIF 305 ! 306 !-- Upper / lower front index for east- and west-facing surfaces 307 ELSEIF ( l == 2 .OR. l == 3 ) THEN 308 IF ( point_index(k,j+1,i) < 0 ) THEN 309 surfaces%npoints = surfaces%npoints + 1 310 point_index(k,j+1,i) = surfaces%npoints - 1 311 ENDIF 312 IF ( point_index(k+1,j+1,i) < 0 ) THEN 313 surfaces%npoints = surfaces%npoints + 1 314 point_index(k+1,j+1,i) = surfaces%npoints - 1 315 ENDIF 316 ENDIF 317 ! 270 318 !-- Upper left / front vertex 271 IF ( point_index(k+1,j,i) < 0 ) THEN319 IF ( point_index(k+1,j,i) < 0 ) THEN 272 320 surfaces%npoints = surfaces%npoints + 1 273 point_index(k+1,j,i) = surfaces%npoints 321 point_index(k+1,j,i) = surfaces%npoints - 1 274 322 ENDIF 323 ENDDO 324 DO m = 1, surf_lsm_v(l)%ns 325 i = surf_lsm_v(l)%i(m) + MERGE( surf_lsm_v(l)%ioff, 0, l == 3 ) 326 j = surf_lsm_v(l)%j(m) + MERGE( surf_lsm_v(l)%joff, 0, l == 1 ) 327 k = surf_lsm_v(l)%k(m) + surf_lsm_v(l)%koff 328 ! 329 !-- Lower left /front vertex 330 IF ( point_index(k,j,i) < 0 ) THEN 331 surfaces%npoints = surfaces%npoints + 1 332 point_index(k,j,i) = surfaces%npoints - 1 333 ENDIF 275 334 ! 276 335 !-- Upper / lower right index for north- and south-facing surfaces 277 IF ( l == 0 . AND. l == 1 ) THEN278 IF ( point_index(k,j,i+1) < 0 ) THEN336 IF ( l == 0 .OR. l == 1 ) THEN 337 IF ( point_index(k,j,i+1) < 0 ) THEN 279 338 surfaces%npoints = surfaces%npoints + 1 280 point_index(k,j,i+1) = surfaces%npoints 339 point_index(k,j,i+1) = surfaces%npoints - 1 281 340 ENDIF 282 IF ( point_index(k+1,j,i+1) < 0 ) THEN341 IF ( point_index(k+1,j,i+1) < 0 ) THEN 283 342 surfaces%npoints = surfaces%npoints + 1 284 point_index(k+1,j,i+1) = surfaces%npoints 343 point_index(k+1,j,i+1) = surfaces%npoints - 1 285 344 ENDIF 286 345 ! 287 346 !-- Upper / lower front index for east- and west-facing surfaces 288 ELSEIF ( l == 2 . AND. l == 3 ) THEN289 IF ( point_index(k,j+1,i) < 0 ) THEN347 ELSEIF ( l == 2 .OR. l == 3 ) THEN 348 IF ( point_index(k,j+1,i) < 0 ) THEN 290 349 surfaces%npoints = surfaces%npoints + 1 291 point_index(k,j+1,i) = surfaces%npoints 350 point_index(k,j+1,i) = surfaces%npoints - 1 292 351 ENDIF 293 IF ( point_index(k+1,j+1,i) < 0 ) THEN352 IF ( point_index(k+1,j+1,i) < 0 ) THEN 294 353 surfaces%npoints = surfaces%npoints + 1 295 point_index(k+1,j+1,i) = surfaces%npoints 296 ENDIF 297 ENDIF 298 299 ENDDO 300 DO m = 1, surf_lsm_v(l)%ns 301 i = surf_lsm_v(l)%i(m) + MERGE( surf_lsm_v(l)%ioff, 0, l == 2 ) 302 j = surf_lsm_v(l)%j(m) + MERGE( surf_lsm_v(l)%joff, 0, l == 0 ) 303 k = surf_lsm_v(l)%k(m) + surf_lsm_v(l)%koff 354 point_index(k+1,j+1,i) = surfaces%npoints - 1 355 ENDIF 356 ENDIF 357 ! 358 !-- Upper left / front vertex 359 IF ( point_index(k+1,j,i) < 0 ) THEN 360 surfaces%npoints = surfaces%npoints + 1 361 point_index(k+1,j,i) = surfaces%npoints - 1 362 ENDIF 363 ENDDO 364 365 DO m = 1, surf_usm_v(l)%ns 366 i = surf_usm_v(l)%i(m) + MERGE( surf_usm_v(l)%ioff, 0, l == 3 ) 367 j = surf_usm_v(l)%j(m) + MERGE( surf_usm_v(l)%joff, 0, l == 1 ) 368 k = surf_usm_v(l)%k(m) + surf_usm_v(l)%koff 304 369 ! 305 370 !-- Lower left /front vertex 306 IF ( point_index(k,j,i) < 0 ) THEN371 IF ( point_index(k,j,i) < 0 ) THEN 307 372 surfaces%npoints = surfaces%npoints + 1 308 point_index(k,j,i) = surfaces%npoints 373 point_index(k,j,i) = surfaces%npoints - 1 309 374 ENDIF 310 375 ! 376 !-- Upper / lower right index for north- and south-facing surfaces 377 IF ( l == 0 .OR. l == 1 ) THEN 378 IF ( point_index(k,j,i+1) < 0 ) THEN 379 surfaces%npoints = surfaces%npoints + 1 380 point_index(k,j,i+1) = surfaces%npoints - 1 381 ENDIF 382 IF ( point_index(k+1,j,i+1) < 0 ) THEN 383 surfaces%npoints = surfaces%npoints + 1 384 point_index(k+1,j,i+1) = surfaces%npoints - 1 385 ENDIF 386 ! 387 !-- Upper / lower front index for east- and west-facing surfaces 388 ELSEIF ( l == 2 .OR. l == 3 ) THEN 389 IF ( point_index(k,j+1,i) < 0 ) THEN 390 surfaces%npoints = surfaces%npoints + 1 391 point_index(k,j+1,i) = surfaces%npoints - 1 392 ENDIF 393 IF ( point_index(k+1,j+1,i) < 0 ) THEN 394 surfaces%npoints = surfaces%npoints + 1 395 point_index(k+1,j+1,i) = surfaces%npoints - 1 396 ENDIF 397 ENDIF 398 ! 311 399 !-- Upper left / front vertex 312 IF ( point_index(k+1,j,i) < 0 ) THEN400 IF ( point_index(k+1,j,i) < 0 ) THEN 313 401 surfaces%npoints = surfaces%npoints + 1 314 point_index(k+1,j,i) = surfaces%npoints 402 point_index(k+1,j,i) = surfaces%npoints - 1 315 403 ENDIF 316 ! 317 !-- Upper / lower right index for north- and south-facing surfaces 318 IF ( l == 0 .AND. l == 1 ) THEN 319 IF ( point_index(k,j,i+1) < 0) THEN 320 surfaces%npoints = surfaces%npoints + 1 321 point_index(k,j,i+1) = surfaces%npoints 322 ENDIF 323 IF ( point_index(k+1,j,i+1) < 0) THEN 324 surfaces%npoints = surfaces%npoints + 1 325 point_index(k+1,j,i+1) = surfaces%npoints 326 ENDIF 327 ! 328 !-- Upper / lower front index for east- and west-facing surfaces 329 ELSEIF ( l == 2 .AND. l == 3 ) THEN 330 IF ( point_index(k,j+1,i) < 0) THEN 331 surfaces%npoints = surfaces%npoints + 1 332 point_index(k,j+1,i) = surfaces%npoints 333 ENDIF 334 IF ( point_index(k+1,j+1,i) < 0) THEN 335 surfaces%npoints = surfaces%npoints + 1 336 point_index(k+1,j+1,i) = surfaces%npoints 337 ENDIF 338 ENDIF 339 ENDDO 340 DO m = 1, surf_usm_v(l)%ns 341 i = surf_usm_v(l)%i(m) + MERGE( surf_usm_v(l)%ioff, 0, l == 2 ) 342 j = surf_usm_v(l)%j(m) + MERGE( surf_usm_v(l)%joff, 0, l == 0 ) 343 k = surf_usm_v(l)%k(m) + surf_usm_v(l)%koff 344 ! 345 !-- Lower left /front vertex 346 IF ( point_index(k,j,i) < 0) THEN 347 surfaces%npoints = surfaces%npoints + 1 348 point_index(k,j,i) = surfaces%npoints 349 ENDIF 350 ! 351 !-- Upper left / front vertex 352 IF ( point_index(k+1,j,i) < 0) THEN 353 surfaces%npoints = surfaces%npoints + 1 354 point_index(k+1,j,i) = surfaces%npoints 355 ENDIF 356 ! 357 !-- Upper / lower right index for north- and south-facing surfaces 358 IF ( l == 0 .AND. l == 1 ) THEN 359 IF ( point_index(k,j,i+1) < 0) THEN 360 surfaces%npoints = surfaces%npoints + 1 361 point_index(k,j,i+1) = surfaces%npoints 362 ENDIF 363 IF ( point_index(k+1,j,i+1) < 0) THEN 364 surfaces%npoints = surfaces%npoints + 1 365 point_index(k+1,j,i+1) = surfaces%npoints 366 ENDIF 367 ! 368 !-- Upper / lower front index for east- and west-facing surfaces 369 ELSEIF ( l == 2 .AND. l == 3 ) THEN 370 IF ( point_index(k,j+1,i) < 0) THEN 371 surfaces%npoints = surfaces%npoints + 1 372 point_index(k,j+1,i) = surfaces%npoints 373 ENDIF 374 IF ( point_index(k+1,j+1,i) < 0) THEN 375 surfaces%npoints = surfaces%npoints + 1 376 point_index(k+1,j+1,i) = surfaces%npoints 377 ENDIF 378 ENDIF 379 ENDDO 404 ENDDO 405 380 406 ENDDO 381 407 ! … … 384 410 !-- of points (vertices) which define the polygons can be larger. 385 411 ALLOCATE( surfaces%points(3,1:surfaces%npoints) ) 386 ALLOCATE( surfaces%polygons(5,1:surfaces%ns) ) 412 ALLOCATE( surfaces%polygons(5,1:surfaces%ns) ) 413 ! 414 !-- Note, PARAVIEW expects consecutively ordered points which can be 415 !-- unambiguously identified. 416 !-- Hence, all PEs know where they should start counting, depending on the 417 !-- number of points on the other PE's with lower MPI rank. 418 #if defined( __parallel ) 419 CALL MPI_ALLGATHER( surfaces%npoints, 1, MPI_INTEGER, & 420 num_points_on_pe, 1, MPI_INTEGER, comm2d, ierr ) 421 #else 422 num_points_on_pe = surfaces%npoints 423 #endif 424 387 425 ! 388 426 !-- After the number of vertices is counted, repeat the loops and define the 389 !-- vertices. 390 !-- Start with the horizontal default surfaces 427 !-- vertices. Start with the horizontal default surfaces 428 !-- First, however, determine the offset where couting of points should be 429 !-- started, which is the sum of points of all PE's with lower MPI rank. 430 i = 0 431 point_index_count = 0 432 DO WHILE ( i < myid .AND. i <= SIZE( num_points_on_pe ) ) 433 point_index_count = point_index_count + num_points_on_pe(i) 434 i = i + 1 435 ENDDO 436 437 391 438 surfaces%npoints = 0 392 439 point_index = -1 393 440 npg = 0 441 394 442 DO l = 0, 1 395 443 DO m = 1, surf_def_h(0)%ns … … 402 450 !-- Check if the vertices that define the surface element are already 403 451 !-- defined, if not, increment the counter. 404 IF ( point_index(k,j,i) < 0 ) THEN452 IF ( point_index(k,j,i) < 0 ) THEN 405 453 surfaces%npoints = surfaces%npoints + 1 406 point_index(k,j,i) = surfaces%npoints 454 point_index(k,j,i) = point_index_count 455 point_index_count = point_index_count + 1 407 456 surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx 408 457 surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy 409 458 surfaces%points(3,surfaces%npoints) = zw(k) 410 459 ENDIF 411 IF ( point_index(k,j +1,i) < 0) THEN460 IF ( point_index(k,j,i+1) < 0 ) THEN 412 461 surfaces%npoints = surfaces%npoints + 1 413 point_index(k,j+1,i) = surfaces%npoints 462 point_index(k,j,i+1) = point_index_count 463 point_index_count = point_index_count + 1 464 surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx 465 surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy 466 surfaces%points(3,surfaces%npoints) = zw(k) 467 ENDIF 468 IF ( point_index(k,j+1,i+1) < 0 ) THEN 469 surfaces%npoints = surfaces%npoints + 1 470 point_index(k,j+1,i+1) = point_index_count 471 point_index_count = point_index_count + 1 472 surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx 473 surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy 474 surfaces%points(3,surfaces%npoints) = zw(k) 475 ENDIF 476 IF ( point_index(k,j+1,i) < 0 ) THEN 477 surfaces%npoints = surfaces%npoints + 1 478 point_index(k,j+1,i) = point_index_count 479 point_index_count = point_index_count + 1 414 480 surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx 415 481 surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy 416 482 surfaces%points(3,surfaces%npoints) = zw(k) 417 483 ENDIF 418 IF ( point_index(k,j,i+1) < 0) THEN419 surfaces%npoints = surfaces%npoints + 1420 point_index(k,j,i+1) = surfaces%npoints421 surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx422 surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy423 surfaces%points(3,surfaces%npoints) = zw(k)424 ENDIF425 IF ( point_index(k,j+1,i+1) < 0) THEN426 surfaces%npoints = surfaces%npoints + 1427 point_index(k,j+1,i+1) = surfaces%npoints428 surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx429 surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy430 surfaces%points(3,surfaces%npoints) = zw(k)431 ENDIF432 484 433 npg = npg + 1485 npg = npg + 1 434 486 surfaces%polygons(1,npg) = 4 435 487 surfaces%polygons(2,npg) = point_index(k,j,i) 436 surfaces%polygons(3,npg) = point_index(k,j +1,i)488 surfaces%polygons(3,npg) = point_index(k,j,i+1) 437 489 surfaces%polygons(4,npg) = point_index(k,j+1,i+1) 438 surfaces%polygons(5,npg) = point_index(k,j ,i+1)490 surfaces%polygons(5,npg) = point_index(k,j+1,i) 439 491 ENDDO 440 492 ENDDO … … 443 495 j = surf_lsm_h%j(m) + surf_lsm_h%joff 444 496 k = surf_lsm_h%k(m) + surf_lsm_h%koff 445 446 IF ( point_index(k,j,i) < 0) THEN 497 IF ( point_index(k,j,i) < 0 ) THEN 447 498 surfaces%npoints = surfaces%npoints + 1 448 point_index(k,j,i) = surfaces%npoints 499 point_index(k,j,i) = point_index_count 500 point_index_count = point_index_count + 1 449 501 surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx 450 502 surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy 451 503 surfaces%points(3,surfaces%npoints) = zw(k) 452 504 ENDIF 453 IF ( point_index(k,j +1,i) < 0) THEN505 IF ( point_index(k,j,i+1) < 0 ) THEN 454 506 surfaces%npoints = surfaces%npoints + 1 455 point_index(k,j+1,i) = surfaces%npoints 507 point_index(k,j,i+1) = point_index_count 508 point_index_count = point_index_count + 1 509 surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx 510 surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy 511 surfaces%points(3,surfaces%npoints) = zw(k) 512 ENDIF 513 IF ( point_index(k,j+1,i+1) < 0 ) THEN 514 surfaces%npoints = surfaces%npoints + 1 515 point_index(k,j+1,i+1) = point_index_count 516 point_index_count = point_index_count + 1 517 surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx 518 surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy 519 surfaces%points(3,surfaces%npoints) = zw(k) 520 ENDIF 521 IF ( point_index(k,j+1,i) < 0 ) THEN 522 surfaces%npoints = surfaces%npoints + 1 523 point_index(k,j+1,i) = point_index_count 524 point_index_count = point_index_count + 1 456 525 surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx 457 526 surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy 458 527 surfaces%points(3,surfaces%npoints) = zw(k) 459 528 ENDIF 460 IF ( point_index(k,j,i+1) < 0) THEN 529 530 npg = npg + 1 531 surfaces%polygons(1,npg) = 4 532 surfaces%polygons(2,npg) = point_index(k,j,i) 533 surfaces%polygons(3,npg) = point_index(k,j,i+1) 534 surfaces%polygons(4,npg) = point_index(k,j+1,i+1) 535 surfaces%polygons(5,npg) = point_index(k,j+1,i) 536 ENDDO 537 538 DO m = 1, surf_usm_h%ns 539 i = surf_usm_h%i(m) + surf_usm_h%ioff 540 j = surf_usm_h%j(m) + surf_usm_h%joff 541 k = surf_usm_h%k(m) + surf_usm_h%koff 542 543 IF ( point_index(k,j,i) < 0 ) THEN 544 surfaces%npoints = surfaces%npoints + 1 545 point_index(k,j,i) = point_index_count 546 point_index_count = point_index_count + 1 547 surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx 548 surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy 549 surfaces%points(3,surfaces%npoints) = zw(k) 550 ENDIF 551 IF ( point_index(k,j,i+1) < 0 ) THEN 461 552 surfaces%npoints = surfaces%npoints + 1 462 point_index(k,j,i+1) = surfaces%npoints 553 point_index(k,j,i+1) = point_index_count 554 point_index_count = point_index_count + 1 463 555 surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx 464 556 surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy 465 557 surfaces%points(3,surfaces%npoints) = zw(k) 466 558 ENDIF 467 IF ( point_index(k,j+1,i+1) < 0 ) THEN559 IF ( point_index(k,j+1,i+1) < 0 ) THEN 468 560 surfaces%npoints = surfaces%npoints + 1 469 point_index(k,j+1,i+1) = surfaces%npoints 561 point_index(k,j+1,i+1) = point_index_count 562 point_index_count = point_index_count + 1 470 563 surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx 471 564 surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy 472 565 surfaces%points(3,surfaces%npoints) = zw(k) 473 566 ENDIF 474 475 npg = npg + 1 476 surfaces%polygons(1,npg) = 4 477 surfaces%polygons(2,npg) = point_index(k,j,i) 478 surfaces%polygons(3,npg) = point_index(k,j+1,i) 479 surfaces%polygons(4,npg) = point_index(k,j+1,i+1) 480 surfaces%polygons(5,npg) = point_index(k,j,i+1) 481 ENDDO 482 DO m = 1, surf_usm_h%ns 483 i = surf_usm_h%i(m) + surf_usm_h%ioff 484 j = surf_usm_h%j(m) + surf_usm_h%joff 485 k = surf_usm_h%k(m) + surf_usm_h%koff 486 487 IF ( point_index(k,j,i) < 0) THEN 488 surfaces%npoints = surfaces%npoints + 1 489 point_index(k,j,i) = surfaces%npoints 490 surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx 491 surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy 492 surfaces%points(3,surfaces%npoints) = zw(k) 493 ENDIF 494 IF ( point_index(k,j+1,i) < 0) THEN 567 IF ( point_index(k,j+1,i) < 0 ) THEN 495 568 surfaces%npoints = surfaces%npoints + 1 496 point_index(k,j+1,i) = surfaces%npoints 569 point_index(k,j+1,i) = point_index_count 570 point_index_count = point_index_count + 1 497 571 surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx 498 572 surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy 499 573 surfaces%points(3,surfaces%npoints) = zw(k) 500 574 ENDIF 501 IF ( point_index(k,j,i+1) < 0) THEN502 surfaces%npoints = surfaces%npoints + 1503 point_index(k,j,i+1) = surfaces%npoints504 surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx505 surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy506 surfaces%points(3,surfaces%npoints) = zw(k)507 ENDIF508 IF ( point_index(k,j+1,i+1) < 0) THEN509 surfaces%npoints = surfaces%npoints + 1510 point_index(k,j+1,i+1) = surfaces%npoints511 surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx512 surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy513 surfaces%points(3,surfaces%npoints) = zw(k)514 ENDIF515 575 516 npg = npg + 1576 npg = npg + 1 517 577 surfaces%polygons(1,npg) = 4 518 578 surfaces%polygons(2,npg) = point_index(k,j,i) 519 surfaces%polygons(3,npg) = point_index(k,j +1,i)579 surfaces%polygons(3,npg) = point_index(k,j,i+1) 520 580 surfaces%polygons(4,npg) = point_index(k,j+1,i+1) 521 surfaces%polygons(5,npg) = point_index(k,j ,i+1)581 surfaces%polygons(5,npg) = point_index(k,j+1,i) 522 582 ENDDO 583 523 584 DO l = 0, 3 524 585 DO m = 1, surf_def_v(l)%ns … … 528 589 !-- identical to the reference j-index outside the grid box. 529 590 !-- Equivalent for east-facing surfaces and i-index. 530 i = surf_def_v(l)%i(m) + MERGE( surf_def_v(l)%ioff, 0, l == 2)531 j = surf_def_v(l)%j(m) + MERGE( surf_def_v(l)%joff, 0, l == 0)591 i = surf_def_v(l)%i(m) + MERGE( surf_def_v(l)%ioff, 0, l == 3 ) 592 j = surf_def_v(l)%j(m) + MERGE( surf_def_v(l)%joff, 0, l == 1 ) 532 593 k = surf_def_v(l)%k(m) + surf_def_v(l)%koff 533 594 ! 534 595 !-- Lower left /front vertex 535 IF ( point_index(k,j,i) < 0 ) THEN596 IF ( point_index(k,j,i) < 0 ) THEN 536 597 surfaces%npoints = surfaces%npoints + 1 537 point_index(k,j,i) = surfaces%npoints 598 point_index(k,j,i) = point_index_count 599 point_index_count = point_index_count + 1 600 surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx 601 surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy 602 surfaces%points(3,surfaces%npoints) = zw(k-1) 603 ENDIF 604 ! 605 !-- Upper / lower right index for north- and south-facing surfaces 606 IF ( l == 0 .OR. l == 1 ) THEN 607 IF ( point_index(k,j,i+1) < 0 ) THEN 608 surfaces%npoints = surfaces%npoints + 1 609 point_index(k,j,i+1) = point_index_count 610 point_index_count = point_index_count + 1 611 surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx 612 surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy 613 surfaces%points(3,surfaces%npoints) = zw(k-1) 614 ENDIF 615 IF ( point_index(k+1,j,i+1) < 0 ) THEN 616 surfaces%npoints = surfaces%npoints + 1 617 point_index(k+1,j,i+1) = point_index_count 618 point_index_count = point_index_count + 1 619 surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx 620 surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy 621 surfaces%points(3,surfaces%npoints) = zw(k) 622 ENDIF 623 ! 624 !-- Upper / lower front index for east- and west-facing surfaces 625 ELSEIF ( l == 2 .OR. l == 3 ) THEN 626 IF ( point_index(k,j+1,i) < 0 ) THEN 627 surfaces%npoints = surfaces%npoints + 1 628 point_index(k,j+1,i) = point_index_count 629 point_index_count = point_index_count + 1 630 surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx 631 surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy 632 surfaces%points(3,surfaces%npoints) = zw(k-1) 633 ENDIF 634 IF ( point_index(k+1,j+1,i) < 0 ) THEN 635 surfaces%npoints = surfaces%npoints + 1 636 point_index(k+1,j+1,i) = point_index_count 637 point_index_count = point_index_count + 1 638 surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx 639 surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy 640 surfaces%points(3,surfaces%npoints) = zw(k) 641 ENDIF 642 ENDIF 643 ! 644 !-- Upper left / front vertex 645 IF ( point_index(k+1,j,i) < 0 ) THEN 646 surfaces%npoints = surfaces%npoints + 1 647 point_index(k+1,j,i) = point_index_count 648 point_index_count = point_index_count + 1 538 649 surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx 539 650 surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy 540 651 surfaces%points(3,surfaces%npoints) = zw(k) 541 ENDIF542 !543 !-- Upper left / front vertex544 IF ( point_index(k+1,j,i) < 0) THEN545 surfaces%npoints = surfaces%npoints + 1546 point_index(k+1,j,i) = surfaces%npoints547 surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx548 surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy549 surfaces%points(3,surfaces%npoints) = zw(k+1)550 652 ENDIF 551 !552 !-- Upper / lower right index for north- and south-facing surfaces553 IF ( l == 0 .AND. l == 1 ) THEN554 IF ( point_index(k,j,i+1) < 0) THEN555 surfaces%npoints = surfaces%npoints + 1556 point_index(k,j,i+1) = surfaces%npoints557 surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx558 surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy559 surfaces%points(3,surfaces%npoints) = zw(k)560 ENDIF561 IF ( point_index(k+1,j,i+1) < 0) THEN562 surfaces%npoints = surfaces%npoints + 1563 point_index(k+1,j,i+1) = surfaces%npoints564 surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx565 surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy566 surfaces%points(3,surfaces%npoints) = zw(k+1)567 ENDIF568 !569 !-- Upper / lower front index for east- and west-facing surfaces570 ELSEIF ( l == 2 .AND. l == 3 ) THEN571 IF ( point_index(k,j+1,i) < 0) THEN572 surfaces%npoints = surfaces%npoints + 1573 point_index(k,j+1,i) = surfaces%npoints574 surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx575 surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy576 surfaces%points(3,surfaces%npoints) = zw(k)577 ENDIF578 IF ( point_index(k+1,j+1,i) < 0) THEN579 surfaces%npoints = surfaces%npoints + 1580 point_index(k+1,j+1,i) = surfaces%npoints581 surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx582 surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy583 surfaces%points(3,surfaces%npoints) = zw(k+1)584 ENDIF585 ENDIF586 653 587 654 npg = npg + 1 588 IF ( l == 0 . AND. l == 1 ) THEN655 IF ( l == 0 .OR. l == 1 ) THEN 589 656 surfaces%polygons(1,npg) = 4 590 657 surfaces%polygons(2,npg) = point_index(k,j,i) … … 601 668 602 669 ENDDO 670 603 671 DO m = 1, surf_lsm_v(l)%ns 604 i = surf_lsm_v(l)%i(m) + MERGE( surf_lsm_v(l)%ioff, 0, l == 2)605 j = surf_lsm_v(l)%j(m) + MERGE( surf_lsm_v(l)%joff, 0, l == 0)672 i = surf_lsm_v(l)%i(m) + MERGE( surf_lsm_v(l)%ioff, 0, l == 3 ) 673 j = surf_lsm_v(l)%j(m) + MERGE( surf_lsm_v(l)%joff, 0, l == 1 ) 606 674 k = surf_lsm_v(l)%k(m) + surf_lsm_v(l)%koff 607 675 ! 608 676 !-- Lower left /front vertex 609 IF ( point_index(k,j,i) < 0 ) THEN677 IF ( point_index(k,j,i) < 0 ) THEN 610 678 surfaces%npoints = surfaces%npoints + 1 611 point_index(k,j,i) = surfaces%npoints 679 point_index(k,j,i) = point_index_count 680 point_index_count = point_index_count + 1 681 surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx 682 surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy 683 surfaces%points(3,surfaces%npoints) = zw(k-1) 684 ENDIF 685 ! 686 !-- Upper / lower right index for north- and south-facing surfaces 687 IF ( l == 0 .OR. l == 1 ) THEN 688 IF ( point_index(k,j,i+1) < 0 ) THEN 689 surfaces%npoints = surfaces%npoints + 1 690 point_index(k,j,i+1) = point_index_count 691 point_index_count = point_index_count + 1 692 surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx 693 surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy 694 surfaces%points(3,surfaces%npoints) = zw(k-1) 695 ENDIF 696 IF ( point_index(k+1,j,i+1) < 0 ) THEN 697 surfaces%npoints = surfaces%npoints + 1 698 point_index(k+1,j,i+1) = point_index_count 699 point_index_count = point_index_count + 1 700 surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx 701 surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy 702 surfaces%points(3,surfaces%npoints) = zw(k) 703 ENDIF 704 ! 705 !-- Upper / lower front index for east- and west-facing surfaces 706 ELSEIF ( l == 2 .OR. l == 3 ) THEN 707 IF ( point_index(k,j+1,i) < 0 ) THEN 708 surfaces%npoints = surfaces%npoints + 1 709 point_index(k,j+1,i) = point_index_count 710 point_index_count = point_index_count + 1 711 surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx 712 surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy 713 surfaces%points(3,surfaces%npoints) = zw(k-1) 714 ENDIF 715 IF ( point_index(k+1,j+1,i) < 0 ) THEN 716 surfaces%npoints = surfaces%npoints + 1 717 point_index(k+1,j+1,i) = point_index_count 718 point_index_count = point_index_count + 1 719 surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx 720 surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy 721 surfaces%points(3,surfaces%npoints) = zw(k) 722 ENDIF 723 ENDIF 724 ! 725 !-- Upper left / front vertex 726 IF ( point_index(k+1,j,i) < 0 ) THEN 727 surfaces%npoints = surfaces%npoints + 1 728 point_index(k+1,j,i) = point_index_count 729 point_index_count = point_index_count + 1 612 730 surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx 613 731 surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy 614 732 surfaces%points(3,surfaces%npoints) = zw(k) 615 ENDIF616 !617 !-- Upper left / front vertex618 IF ( point_index(k+1,j,i) < 0) THEN619 surfaces%npoints = surfaces%npoints + 1620 point_index(k+1,j,i) = surfaces%npoints621 surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx622 surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy623 surfaces%points(3,surfaces%npoints) = zw(k+1)624 ENDIF625 !626 !-- Upper / lower right index for north- and south-facing surfaces627 IF ( l == 0 .AND. l == 1 ) THEN628 IF ( point_index(k,j,i+1) < 0) THEN629 surfaces%npoints = surfaces%npoints + 1630 point_index(k,j,i+1) = surfaces%npoints631 surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx632 surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy633 surfaces%points(3,surfaces%npoints) = zw(k)634 ENDIF635 IF ( point_index(k+1,j,i+1) < 0) THEN636 surfaces%npoints = surfaces%npoints + 1637 point_index(k+1,j,i+1) = surfaces%npoints638 surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx639 surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy640 surfaces%points(3,surfaces%npoints) = zw(k+1)641 ENDIF642 !643 !-- Upper / lower front index for east- and west-facing surfaces644 ELSEIF ( l == 2 .AND. l == 3 ) THEN645 IF ( point_index(k,j+1,i) < 0) THEN646 surfaces%npoints = surfaces%npoints + 1647 point_index(k,j+1,i) = surfaces%npoints648 surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx649 surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy650 surfaces%points(3,surfaces%npoints) = zw(k)651 ENDIF652 IF ( point_index(k+1,j+1,i) < 0) THEN653 surfaces%npoints = surfaces%npoints + 1654 point_index(k+1,j+1,i) = surfaces%npoints655 surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx656 surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy657 surfaces%points(3,surfaces%npoints) = zw(k+1)658 ENDIF659 733 ENDIF 660 734 661 735 npg = npg + 1 662 IF ( l == 0 . AND. l == 1 ) THEN736 IF ( l == 0 .OR. l == 1 ) THEN 663 737 surfaces%polygons(1,npg) = 4 664 738 surfaces%polygons(2,npg) = point_index(k,j,i) … … 675 749 ENDDO 676 750 DO m = 1, surf_usm_v(l)%ns 677 i = surf_usm_v(l)%i(m) + MERGE( surf_usm_v(l)%ioff, 0, l == 2)678 j = surf_usm_v(l)%j(m) + MERGE( surf_usm_v(l)%joff, 0, l == 0)751 i = surf_usm_v(l)%i(m) + MERGE( surf_usm_v(l)%ioff, 0, l == 3 ) 752 j = surf_usm_v(l)%j(m) + MERGE( surf_usm_v(l)%joff, 0, l == 1 ) 679 753 k = surf_usm_v(l)%k(m) + surf_usm_v(l)%koff 680 754 ! 681 755 !-- Lower left /front vertex 682 IF ( point_index(k,j,i) < 0 ) THEN756 IF ( point_index(k,j,i) < 0 ) THEN 683 757 surfaces%npoints = surfaces%npoints + 1 684 point_index(k,j,i) = surfaces%npoints 758 point_index(k,j,i) = point_index_count 759 point_index_count = point_index_count + 1 760 surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx 761 surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy 762 surfaces%points(3,surfaces%npoints) = zw(k-1) 763 ENDIF 764 ! 765 !-- Upper / lower right index for north- and south-facing surfaces 766 IF ( l == 0 .OR. l == 1 ) THEN 767 IF ( point_index(k,j,i+1) < 0 ) THEN 768 surfaces%npoints = surfaces%npoints + 1 769 point_index(k,j,i+1) = point_index_count 770 point_index_count = point_index_count + 1 771 surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx 772 surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy 773 surfaces%points(3,surfaces%npoints) = zw(k-1) 774 ENDIF 775 IF ( point_index(k+1,j,i+1) < 0 ) THEN 776 surfaces%npoints = surfaces%npoints + 1 777 point_index(k+1,j,i+1) = point_index_count 778 point_index_count = point_index_count + 1 779 surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx 780 surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy 781 surfaces%points(3,surfaces%npoints) = zw(k) 782 ENDIF 783 ! 784 !-- Upper / lower front index for east- and west-facing surfaces 785 ELSEIF ( l == 2 .OR. l == 3 ) THEN 786 IF ( point_index(k,j+1,i) < 0 ) THEN 787 surfaces%npoints = surfaces%npoints + 1 788 point_index(k,j+1,i) = point_index_count 789 point_index_count = point_index_count + 1 790 surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx 791 surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy 792 surfaces%points(3,surfaces%npoints) = zw(k-1) 793 ENDIF 794 IF ( point_index(k+1,j+1,i) < 0 ) THEN 795 surfaces%npoints = surfaces%npoints + 1 796 point_index(k+1,j+1,i) = point_index_count 797 point_index_count = point_index_count + 1 798 surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx 799 surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy 800 surfaces%points(3,surfaces%npoints) = zw(k) 801 ENDIF 802 ENDIF 803 ! 804 !-- Upper left / front vertex 805 IF ( point_index(k+1,j,i) < 0 ) THEN 806 surfaces%npoints = surfaces%npoints + 1 807 point_index(k+1,j,i) = point_index_count 808 point_index_count = point_index_count + 1 685 809 surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx 686 810 surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy 687 811 surfaces%points(3,surfaces%npoints) = zw(k) 688 ENDIF689 !690 !-- Upper left / front vertex691 IF ( point_index(k+1,j,i) < 0) THEN692 surfaces%npoints = surfaces%npoints + 1693 point_index(k+1,j,i) = surfaces%npoints694 surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx695 surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy696 surfaces%points(3,surfaces%npoints) = zw(k+1)697 ENDIF698 !699 !-- Upper / lower right index for north- and south-facing surfaces700 IF ( l == 0 .AND. l == 1 ) THEN701 IF ( point_index(k,j,i+1) < 0) THEN702 surfaces%npoints = surfaces%npoints + 1703 point_index(k,j,i+1) = surfaces%npoints704 surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx705 surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy706 surfaces%points(3,surfaces%npoints) = zw(k)707 ENDIF708 IF ( point_index(k+1,j,i+1) < 0) THEN709 surfaces%npoints = surfaces%npoints + 1710 point_index(k+1,j,i+1) = surfaces%npoints711 surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx712 surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy713 surfaces%points(3,surfaces%npoints) = zw(k+1)714 ENDIF715 !716 !-- Upper / lower front index for east- and west-facing surfaces717 ELSEIF ( l == 2 .AND. l == 3 ) THEN718 IF ( point_index(k,j+1,i) < 0) THEN719 surfaces%npoints = surfaces%npoints + 1720 point_index(k,j+1,i) = surfaces%npoints721 surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx722 surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy723 surfaces%points(3,surfaces%npoints) = zw(k)724 ENDIF725 IF ( point_index(k+1,j+1,i) < 0) THEN726 surfaces%npoints = surfaces%npoints + 1727 point_index(k+1,j+1,i) = surfaces%npoints728 surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx729 surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy730 surfaces%points(3,surfaces%npoints) = zw(k+1)731 ENDIF732 812 ENDIF 733 813 734 814 npg = npg + 1 735 IF ( l == 0 . AND. l == 1 ) THEN815 IF ( l == 0 .OR. l == 1 ) THEN 736 816 surfaces%polygons(1,npg) = 4 737 817 surfaces%polygons(2,npg) = point_index(k,j,i) … … 747 827 ENDIF 748 828 ENDDO 829 749 830 ENDDO 750 831 ! … … 752 833 DEALLOCATE ( point_index ) 753 834 ! 754 !-- Now, vertices as well as polygon data can be written to NetCDF files 755 !-- ... 835 !-- Sum-up total number of surface elements and vertices on domain. This 836 !-- will be needed for post-processing. 837 surfaces%npoints_total = 0 838 #if defined( __parallel ) 839 CALL MPI_ALLREDUCE( surfaces%npoints, surfaces%npoints_total, 1, & 840 MPI_INTEGER, MPI_SUM, comm2d, ierr ) 841 CALL MPI_ALLREDUCE( surfaces%ns, surfaces%ns_total, 1, & 842 MPI_INTEGER, MPI_SUM, comm2d, ierr ) 843 #else 844 surfaces%npoints_total = surfaces%npoints 845 surfaces%ns_total = surfaces%ns 846 #endif 847 756 848 END SUBROUTINE surface_output_init 757 849 … … 795 887 IF ( i == io_group ) THEN 796 888 WRITE ( 25+av ) surfaces%npoints 889 WRITE ( 25+av ) surfaces%npoints_total 797 890 WRITE ( 25+av ) surfaces%ns 891 WRITE ( 25+av ) surfaces%ns_total 798 892 WRITE ( 25+av ) surfaces%points 799 893 WRITE ( 25+av ) surfaces%polygons … … 805 899 ENDDO 806 900 ENDIF 807 !808 !-- Write time coordinate809 DO i = 0, io_blocks-1810 IF ( i == io_group ) THEN811 trimvar = 'time'812 WRITE ( 25+av ) trimvar813 WRITE ( 25+av ) simulated_time814 ENDIF815 #if defined( __parallel )816 CALL MPI_BARRIER( comm2d, ierr )817 #endif818 ENDDO819 820 901 821 902 n_out = 0 … … 835 916 IF ( av == 0 ) THEN 836 917 CALL surface_output_collect( surf_def_h(0)%us, & 837 surf_def_h(1)%us,&838 surf_lsm_h%us,&839 surf_usm_h%us,&840 surf_def_v(0)%us,&841 surf_lsm_v(0)%us,&842 surf_usm_v(0)%us,&843 surf_def_v(1)%us,&844 surf_lsm_v(1)%us,&845 surf_usm_v(1)%us,&846 surf_def_v(2)%us,&847 surf_lsm_v(2)%us,&848 surf_usm_v(2)%us,&849 surf_def_v(3)%us,&850 surf_lsm_v(3)%us,&851 surf_usm_v(3)%us )918 surf_def_h(1)%us, & 919 surf_lsm_h%us, & 920 surf_usm_h%us, & 921 surf_def_v(0)%us, & 922 surf_lsm_v(0)%us, & 923 surf_usm_v(0)%us, & 924 surf_def_v(1)%us, & 925 surf_lsm_v(1)%us, & 926 surf_usm_v(1)%us, & 927 surf_def_v(2)%us, & 928 surf_lsm_v(2)%us, & 929 surf_usm_v(2)%us, & 930 surf_def_v(3)%us, & 931 surf_lsm_v(3)%us, & 932 surf_usm_v(3)%us ) 852 933 ELSE 853 934 ! … … 1945 2026 DO i = 0, io_blocks-1 1946 2027 IF ( i == io_group ) THEN 1947 WRITE ( 25+av ) trimvar 2028 WRITE ( 25+av ) LEN_TRIM( 'time' ) 2029 WRITE ( 25+av ) 'time' 2030 WRITE ( 25+av ) simulated_time 2031 WRITE ( 25+av ) LEN_TRIM( trimvar ) 2032 WRITE ( 25+av ) TRIM( trimvar ) 1948 2033 WRITE ( 25+av ) surfaces%var_out 1949 2034 ENDIF … … 1956 2041 1957 2042 ! 1958 !-- 1959 2043 !-- If averaged output was written to NetCDF file, set the counter to zero 2044 IF ( av == 1 ) average_count_surf = 0 1960 2045 1961 2046 END SUBROUTINE surface_output … … 2679 2764 !> Sum-up the surface data for average output variables. 2680 2765 !------------------------------------------------------------------------------! 2681 SUBROUTINE surface_output_sum_up( var_def_h0, var_def_h1, & 2682 var_lsm_h, var_usm_h, & 2683 var_def_v0, var_def_v1, var_def_v2, & 2684 var_def_v3, & 2685 var_lsm_v0, var_lsm_v1, var_lsm_v2, & 2686 var_lsm_v3, & 2687 var_usm_v0, var_usm_v1, var_usm_v2, & 2688 var_usm_v3, n_out ) 2766 SUBROUTINE surface_output_sum_up( var_def_h0, var_def_h1, & 2767 var_lsm_h, var_usm_h, & 2768 var_def_v0, var_lsm_v0, var_usm_v0, & 2769 var_def_v1, var_lsm_v1, var_usm_v1, & 2770 var_def_v2, var_lsm_v2, var_usm_v2, & 2771 var_def_v3, var_lsm_v3, var_usm_v3, n_out ) 2689 2772 2690 2773 IMPLICIT NONE … … 2886 2969 !------------------------------------------------------------------------------! 2887 2970 SUBROUTINE surface_output_collect( var_def_h0, var_def_h1, & 2888 var_lsm_h, var_usm_h, & 2889 var_def_v0, var_def_v1, var_def_v2, & 2890 var_def_v3, & 2891 var_lsm_v0, var_lsm_v1, var_lsm_v2, & 2892 var_lsm_v3, & 2893 var_usm_v0, var_usm_v1, var_usm_v2, & 2894 var_usm_v3 ) 2971 var_lsm_h, var_usm_h, & 2972 var_def_v0, var_lsm_v0, var_usm_v0, & 2973 var_def_v1, var_lsm_v1, var_usm_v1, & 2974 var_def_v2, var_lsm_v2, var_usm_v2, & 2975 var_def_v3, var_lsm_v3, var_usm_v3 ) 2895 2976 2896 2977 IMPLICIT NONE … … 3246 3327 3247 3328 END SUBROUTINE surface_output_check_parameters 3329 3330 3331 !------------------------------------------------------------------------------! 3332 ! Description: 3333 ! ------------ 3334 !> Last action. 3335 !------------------------------------------------------------------------------! 3336 SUBROUTINE surface_output_last_action( av ) 3337 3338 USE control_parameters, & 3339 ONLY: io_blocks, io_group 3340 3341 USE pegrid, & 3342 ONLY: comm2d, ierr 3343 3344 IMPLICIT NONE 3345 3346 INTEGER(iwp) :: av !< id indicating average or non-average data output 3347 INTEGER(iwp) :: i !< loop index 3348 3349 ! 3350 !-- Return, if nothing to output 3351 IF ( dosurf_no(av) == 0 ) RETURN 3352 ! 3353 !-- Open files 3354 CALL check_open( 25+av ) 3355 ! 3356 !-- Write time coordinate 3357 DO i = 0, io_blocks-1 3358 IF ( i == io_group ) THEN 3359 WRITE ( 25+av ) LEN_TRIM( 'END' ) 3360 WRITE ( 25+av ) 'END' 3361 ENDIF 3362 #if defined( __parallel ) 3363 CALL MPI_BARRIER( comm2d, ierr ) 3364 #endif 3365 ENDDO 3366 3367 END SUBROUTINE surface_output_last_action 3248 3368 3249 3369 !--Private Subroutines (Only called from tinside his module)
Note: See TracChangeset
for help on using the changeset viewer.
