- Timestamp:
- Nov 4, 2020 1:12:46 PM (4 years ago)
- Location:
- palm/trunk/SOURCE
- Files:
-
- 2 edited
-
pmc_handle_communicator_mod.f90 (modified) (8 diffs)
-
pmc_interface_mod.f90 (modified) (9 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/pmc_handle_communicator_mod.f90
r4650 r4771 25 25 ! ----------------- 26 26 ! $Id$ 27 ! Canopy-restricted anterpolation introduced. New namelist parameter anterpolation_starting_height 28 ! introduced for controlling canopy-restricted anterpolation. 29 ! 30 ! 4650 2020-08-25 14:35:50Z raasch 27 31 ! bugfix for r4649 28 32 ! … … 138 142 !--------------------------------------------------------------------------------------------------! 139 143 SUBROUTINE pmc_init_model( comm, nesting_datatransfer_mode, nesting_mode, & 140 anterpolation_buffer_width, pmc_status )144 anterpolation_buffer_width, anterpolation_starting_height, pmc_status ) 141 145 142 146 USE control_parameters, & … … 155 159 INTEGER, INTENT(INOUT) :: pmc_status !< 156 160 161 REAL(wp), INTENT(INOUT) :: anterpolation_starting_height !< steering parameter for canopy restricted anterpolation 162 157 163 INTEGER :: childcount !< 158 164 INTEGER :: i !< … … 180 186 181 187 CALL read_coupling_layout( nesting_datatransfer_mode, nesting_mode, & 182 anterpolation_buffer_width, pmc_status ) 188 anterpolation_buffer_width, anterpolation_starting_height, & 189 pmc_status ) 183 190 184 191 IF ( pmc_status /= pmc_no_namelist_found .AND. & … … 252 259 MPI_COMM_WORLD, istat ) 253 260 CALL MPI_BCAST( anterpolation_buffer_width, 1, MPI_INT, 0, MPI_COMM_WORLD, istat ) 261 CALL MPI_BCAST( anterpolation_starting_height, 1, MPI_REAL, 0, MPI_COMM_WORLD, istat ) 254 262 ! 255 263 !-- Assign global MPI processes to individual models by setting the couple id … … 422 430 !--------------------------------------------------------------------------------------------------! 423 431 SUBROUTINE read_coupling_layout( nesting_datatransfer_mode, nesting_mode, & 424 anterpolation_buffer_width, pmc_status ) 432 anterpolation_buffer_width, anterpolation_starting_height, & 433 pmc_status ) 425 434 426 435 IMPLICIT NONE … … 431 440 INTEGER, INTENT(INOUT) :: anterpolation_buffer_width !< Boundary buffer width for anterpolation 432 441 INTEGER(iwp), INTENT(INOUT) :: pmc_status !< 442 443 REAL(wp), INTENT(INOUT) :: anterpolation_starting_height !< steering parameter for canopy restricted anterpolation 433 444 434 445 INTEGER(iwp) :: bad_llcorner !< … … 441 452 nesting_datatransfer_mode, & 442 453 nesting_mode, & 443 anterpolation_buffer_width 454 anterpolation_buffer_width, & 455 anterpolation_starting_height 444 456 445 457 ! -
palm/trunk/SOURCE/pmc_interface_mod.f90
r4745 r4771 25 25 ! ----------------- 26 26 ! $Id$ 27 ! Canopy-restricted anterpolation introduced. New namelist parameter anterpolation_starting_height 28 ! introduced for controlling canopy-restricted anterpolation. 29 ! 30 ! 4745 2020-10-15 16:37:13Z suehring 27 31 ! Adjustement of face-area calculation to 3D-topographies 28 32 ! … … 453 457 INTEGER(iwp), SAVE :: anterpolation_buffer_width = 2 !< Boundary buffer width for anterpolation 454 458 459 REAL(wp), SAVE :: anterpolation_starting_height = 9999999.9_wp !< steering parameter for canopy restricted anterpolation 460 455 461 LOGICAL, SAVE :: nested_run = .FALSE. !< general switch 456 462 LOGICAL, SAVE :: rans_mode_parent = .FALSE. !< mode of parent model (.F. - LES mode, .T. - RANS mode) … … 463 469 ! 464 470 !-- Children's parent-grid arrays 471 INTEGER(iwp), SAVE, DIMENSION(:,:), ALLOCATABLE :: kpb_anterp !< Lower limit of kp in anterpolation 465 472 INTEGER(iwp), SAVE, DIMENSION(5), PUBLIC :: parent_bound !< subdomain index bounds for children's parent-grid arrays 466 473 … … 702 709 703 710 CALL pmc_init_model( world_comm, nesting_datatransfer_mode, nesting_mode, & 704 anterpolation_buffer_width, pmc_status )711 anterpolation_buffer_width, anterpolation_starting_height, pmc_status ) 705 712 706 713 IF ( pmc_status == pmc_no_namelist_found ) THEN … … 1507 1514 !-- Compute surface areas of the nest-boundary faces 1508 1515 CALL pmci_compute_face_areas 1509 1516 ! 1517 !-- Compute anterpolation lower kp-index bounds if necessary 1518 CALL pmci_compute_kpb_anterp 1510 1519 ENDIF 1511 1520 … … 2438 2447 2439 2448 END SUBROUTINE pmci_compute_face_areas 2449 2450 2451 !--------------------------------------------------------------------------------------------------! 2452 ! Description: 2453 ! ------------ 2454 !> Define the anterpolation starting kp-indices for all jp,ip based on the obstacle-canopy 2455 !> topograhy for canopy-restricted anterpolation. Note that this is based on the child terrain 2456 !> topography information since it is difficult to access the parent topography information 2457 !> from the child. This means that these topographies are assumed to be close to each other. 2458 !--------------------------------------------------------------------------------------------------! 2459 SUBROUTINE pmci_compute_kpb_anterp 2460 2461 INTEGER(iwp) :: ic !< Child-grid index in the x-direction 2462 INTEGER(iwp) :: ip !< Parent-grid index in the x-direction 2463 INTEGER(iwp) :: jc !< Child-grid index in the y-direction 2464 INTEGER(iwp) :: jp !< Parent-grid index in the y-direction 2465 INTEGER(iwp) :: terrain_surf_k_child !< Local terrain height k index in the child grid 2466 INTEGER(iwp) :: terrain_surf_k_parent !< Local terrain height k index in the parent grid 2467 2468 2469 ALLOCATE( kpb_anterp(jps:jpn,ipl:ipr) ) 2470 kpb_anterp(jps:jpn,ipl:ipr) = 0 2471 2472 IF ( nesting_mode /= 'one-way' .AND. topography /= 'flat' ) THEN 2473 ! 2474 !-- Check if the value was given by user or not 2475 IF ( anterpolation_starting_height > 9000000.0_wp ) THEN 2476 ! 2477 !-- No value was given by user, compute a default value first. 2478 CALL default_anterpolation_starting_height( 99 ) 2479 ENDIF 2480 2481 DO ip = ipl, ipr 2482 DO jp = jps, jpn 2483 terrain_surf_k_parent = 0 2484 DO ic = iflo(ip), ifuo(ip) 2485 DO jc = jflo(jp), jfuo(jp) 2486 terrain_surf_k_child = MINLOC( & 2487 MERGE( 1, 0, BTEST( wall_flags_total_0(:,jc,ic), 5 ) ), DIM = 1 ) - 1 2488 terrain_surf_k_parent = MAX( terrain_surf_k_child, terrain_surf_k_parent ) 2489 ENDDO 2490 ENDDO 2491 terrain_surf_k_parent = NINT( terrain_surf_k_parent / REAL( kgsr, wp ) ) 2492 kpb_anterp(jp,ip) = terrain_surf_k_parent + NINT( anterpolation_starting_height / pg%dz ) 2493 ENDDO 2494 ENDDO 2495 2496 ENDIF 2497 2498 END SUBROUTINE pmci_compute_kpb_anterp 2499 2500 2501 !--------------------------------------------------------------------------------------------------! 2502 ! Description: 2503 ! ------------ 2504 !> Compute a default value for the anterpolation_starting_height for the canopy-restricted 2505 !> anterpolation. The default value is based on a given percentile of the child building height 2506 !> distribution. This operation is based on the child-grid information. 2507 !--------------------------------------------------------------------------------------------------! 2508 SUBROUTINE default_anterpolation_starting_height( percentile_level ) 2509 2510 INTEGER(iwp), INTENT(IN) :: percentile_level !< Selected percentile level (1,...,99) 2511 2512 INTEGER(iwp) :: ic !< Child-grid index in the x-direction 2513 INTEGER(iwp) :: ierr !< MPI-error code 2514 INTEGER(iwp) :: im !< Grid index in the x-direction of the currently 2515 !< globally highest building node 2516 INTEGER(iwp) :: jc !< Child-grid index in the y-direction 2517 INTEGER(iwp), DIMENSION(2) :: ji_maxloc !< jc and ic indices of the local maximum height 2518 INTEGER(iwp) :: jm !< Grid index in the y-direction of the currently 2519 !< globally highest building node 2520 INTEGER(iwp) :: num_building_nodes !< Global number of jc,ic-nodes under buildings 2521 INTEGER(iwp) :: num_building_nodes_l !< Local number of jc,ic-nodes under buildings 2522 INTEGER(iwp) :: num_highest_nodes !< The number of highest building nodes above the 2523 !< percentile value 2524 INTEGER(iwp) :: n !< Running index for the sorting loop 2525 INTEGER(iwp) :: pid_highest !< Process id of the process where the currently 2526 !< highest building node was found. 2527 INTEGER(iwp) :: terrain_surf_k !< Local terrain height k index in the child grid 2528 2529 REAL(wp) :: global_currently_highest !< Globally highest building node currently in the 2530 !< search loop 2531 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: building_height !< Temporary array for buiding heights 2532 REAL(wp), DIMENSION(:), ALLOCATABLE :: max_height !< Maximum buiding heights of each 2533 !< process 2534 REAL(wp), DIMENSION(:), ALLOCATABLE :: max_height_l !< max_height_l( myid ) is the local 2535 !< maximum buiding height of the 2536 !< current process while the other 2537 !< elements are zeroes 2538 2539 ALLOCATE( building_height(nys:nyn,nxl:nxr) ) 2540 ! 2541 !-- Find the building heights and temporarily store them in building_height(nys:nyn,nxl:nxr). 2542 !-- Also find the number of nodes occupied by buildings 2543 num_building_nodes_l = 0 2544 num_building_nodes = 0 2545 DO ic = nxl, nxr 2546 DO jc = nys, nyn 2547 terrain_surf_k = MINLOC( MERGE( 1, 0, BTEST( wall_flags_total_0(:,jc,ic), 5 ) ), & 2548 DIM = 1 ) - 1 2549 IF ( topo_top_ind(jc,ic,5) > terrain_surf_k ) THEN 2550 building_height(jc,ic) = dz(1) * ( topo_top_ind(jc,ic,5) - terrain_surf_k + 1 ) 2551 num_building_nodes_l = num_building_nodes_l + 1 2552 ENDIF 2553 ENDDO 2554 ENDDO 2555 2556 CALL MPI_ALLREDUCE( num_building_nodes_l, num_building_nodes, 1, MPI_INT, MPI_SUM, comm2d, ierr ) 2557 2558 num_highest_nodes = MAX( 1, NINT( num_building_nodes * ( 100 - percentile_level ) / 100.0_wp ) ) 2559 ! 2560 !-- Find global percentile value 2561 ALLOCATE( max_height_l(0:numprocs-1) ) 2562 ALLOCATE( max_height(0:numprocs-1) ) 2563 max_height_l = 0.0_wp 2564 ! 2565 !-- Search loop for the num_highest_nodes highest nodes globally. 2566 DO n = 1, num_highest_nodes 2567 max_height = 0.0_wp 2568 max_height_l(myid) = MAXVAL( building_height ) 2569 CALL MPI_ALLREDUCE( max_height_l, max_height, numprocs, MPI_REAL, MPI_SUM, comm2d, ierr ) 2570 global_currently_highest = MAXVAL( max_height ) 2571 pid_highest = MAXLOC( max_height, DIM = 1 ) - 1 2572 ! 2573 !-- Set the globally currently highest height to zero in order to find the next 2574 !-- highest node in the next loop cycle. 2575 IF ( pid_highest == myid ) THEN 2576 ji_maxloc = MAXLOC( building_height ) 2577 jm = ji_maxloc(1) - 1 + nys 2578 im = ji_maxloc(2) - 1 + nxl 2579 building_height(jm,im) = 0.0_wp 2580 ENDIF 2581 ENDDO 2582 ! 2583 !-- Now after the search-loop, global_currently_highest is the desired percentile 2584 !-- value. Add two parent-grid dz in order to create appropriate clearance above the 2585 !-- roofs in cases of evenly distributed building height such as simple 2586 !-- cuboid-array test cases. 2587 anterpolation_starting_height = global_currently_highest + 2.0_wp * pg%dz 2588 2589 WRITE( 9, * ) 2590 WRITE( 9, "('A default value is defined and set for anterpolation_starting_height = ', & 2591 f4.1, ' m.' )" ) anterpolation_starting_height 2592 WRITE( 9, * ) 2593 2594 DEALLOCATE( max_height_l ) 2595 DEALLOCATE( max_height ) 2596 DEALLOCATE( building_height ) 2597 2598 END SUBROUTINE default_anterpolation_starting_height 2599 2440 2600 #endif 2441 2601 … … 4815 4975 INTEGER(iwp) :: kc !< Running index z-direction - child grid 4816 4976 INTEGER(iwp) :: kp !< Running index z-direction - parent grid 4817 INTEGER(iwp) :: kpb_anterp = 0 !< Bottom boundary index for anterpolation along z4818 4977 INTEGER(iwp) :: kpt_anterp !< Top boundary index for anterpolation along z 4819 4978 INTEGER(iwp) :: var_flag !< bit number used to flag topography on respective grid … … 4834 4993 jps_anterp = jps 4835 4994 jpn_anterp = jpn 4836 kpb_anterp = 0 4995 ! 4996 !-- kpb_anterp is a function of jp and ip and it is set in the initialization phase in 4997 !-- pmci_compute_kpb_anterp. 4837 4998 kpt_anterp = kct - 1 - anterpolation_buffer_width 4838 4999 … … 4862 5023 DO jp = jps_anterp, jpn_anterp 4863 5024 ! 4864 !-- For simplicity anterpolate within buildings and under elevated terrain too 4865 DO kp = kpb_anterp, kpt_anterp 5025 !-- If the user has set anterpolation_starting_height less than the canopy height, the 5026 !-- anterpolation is made also within buildings for simplicity, and even under elevated 5027 !-- terrain if anterpolation_starting_height is set smaller than terrain height. 5028 DO kp = kpb_anterp(jp,ip), kpt_anterp 4866 5029 cellsum = 0.0_wp 4867 5030 DO ic = ifl(ip), ifu(ip)
Note: See TracChangeset
for help on using the changeset viewer.
