Changeset 2232 for palm/trunk/SOURCE/pres.f90
- Timestamp:
- May 30, 2017 5:47:52 PM (7 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/pres.f90 (modified) (26 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/pres.f90
r2119 r2232 20 20 ! Current revisions: 21 21 ! ------------------ 22 ! 22 ! Adjustments to new topography and surface concept 23 23 ! 24 24 ! Former revisions: … … 156 156 USE indices, & 157 157 ONLY: nbgp, ngp_2dh_outer, nx, nxl, nxlg, nxl_mg, nxr, nxrg, nxr_mg, & 158 ny, nys, nysg, nys_mg, nyn, nyng, nyn_mg, nzb, nzb_s_inner, & 159 nzb_u_inner, nzb_v_inner, nzb_w_inner, nzt, nzt_mg, & 160 rflags_s_inner 158 ny, nys, nysg, nys_mg, nyn, nyng, nyn_mg, nzb, nzt, nzt_mg, & 159 wall_flags_0 161 160 162 161 USE kinds … … 176 175 weight_substep 177 176 177 USE surface_mod, & 178 ONLY : bc_h 179 178 180 IMPLICIT NONE 179 181 … … 181 183 INTEGER(iwp) :: j !< 182 184 INTEGER(iwp) :: k !< 185 INTEGER(iwp) :: m !< 183 186 184 187 REAL(wp) :: ddt_3d !< … … 269 272 ! 270 273 !-- Sum up the volume flow through the south/north boundary 271 DO k = nzb_u_inner(j,i)+1, nzt 272 volume_flow_l(1) = volume_flow_l(1) + u(k,j,i) * dzw(k) 274 DO k = nzb+1, nzt 275 volume_flow_l(1) = volume_flow_l(1) + u(k,j,i) * dzw(k) & 276 * MERGE( 1.0_wp, 0.0_wp, & 277 BTEST( wall_flags_0(k,j,i), 1 ) & 278 ) 273 279 ENDDO 274 280 ENDDO … … 285 291 286 292 DO j = nysg, nyng 287 DO k = nzb_u_inner(j,i)+1, nzt 288 u(k,j,i) = u(k,j,i) + volume_flow_offset(1) 293 DO k = nzb+1, nzt 294 u(k,j,i) = u(k,j,i) + volume_flow_offset(1) & 295 * MERGE( 1.0_wp, 0.0_wp, & 296 BTEST( wall_flags_0(k,j,i), 1 ) & 297 ) 289 298 ENDDO 290 299 ENDDO … … 308 317 ! 309 318 !-- Sum up the volume flow through the south/north boundary 310 DO k = nzb_v_inner(j,i)+1, nzt 311 volume_flow_l(2) = volume_flow_l(2) + v(k,j,i) * dzw(k) 319 DO k = nzb+1, nzt 320 volume_flow_l(2) = volume_flow_l(2) + v(k,j,i) * dzw(k) & 321 * MERGE( 1.0_wp, 0.0_wp, & 322 BTEST( wall_flags_0(k,j,i), 2 ) & 323 ) 312 324 ENDDO 313 325 ENDDO … … 324 336 325 337 DO i = nxlg, nxrg 326 DO k = nzb_v_inner(j,i)+1, nzt 327 v(k,j,i) = v(k,j,i) + volume_flow_offset(2) 338 DO k = nzb+1, nzt 339 v(k,j,i) = v(k,j,i) + volume_flow_offset(2) & 340 * MERGE( 1.0_wp, 0.0_wp, & 341 BTEST( wall_flags_0(k,j,i), 2 ) & 342 ) 328 343 ENDDO 329 344 ENDDO … … 350 365 DO i = nxl, nxr 351 366 DO j = nys, nyn 352 DO k = nzb_w_inner(j,i)+1, nzt 353 w_l_l(k) = w_l_l(k) + w(k,j,i) 367 DO k = nzb+1, nzt 368 w_l_l(k) = w_l_l(k) + w(k,j,i) & 369 * MERGE( 1.0_wp, 0.0_wp, & 370 BTEST( wall_flags_0(k,j,i), 3 ) & 371 ) 354 372 ENDDO 355 373 ENDDO … … 367 385 DO i = nxlg, nxrg 368 386 DO j = nysg, nyng 369 DO k = nzb_w_inner(j,i)+1, nzt 370 w(k,j,i) = w(k,j,i) - w_l(k) 387 DO k = nzb+1, nzt 388 w(k,j,i) = w(k,j,i) - w_l(k) & 389 * MERGE( 1.0_wp, 0.0_wp, & 390 BTEST( wall_flags_0(k,j,i), 3 ) & 391 ) 371 392 ENDDO 372 393 ENDDO … … 379 400 380 401 IF ( psolver(1:9) == 'multigrid' ) THEN 381 !$OMP PARALLEL DO SCHEDULE( STATIC ) 402 !$OMP PARALLEL DO SCHEDULE( STATIC ) PRIVATE (i,j,k) 382 403 DO i = nxl-1, nxr+1 383 404 DO j = nys-1, nyn+1 … … 388 409 ENDDO 389 410 ELSE 390 !$OMP PARALLEL DO SCHEDULE( STATIC ) 411 !$OMP PARALLEL DO SCHEDULE( STATIC ) PRIVATE (i,j,k) 391 412 DO i = nxl, nxr 392 413 DO j = nys, nyn … … 406 427 DO i = nxl, nxr 407 428 DO j = nys, nyn 408 DO k = nzb _s_inner(j,i)+1, nzt429 DO k = nzb+1, nzt 409 430 d(k,j,i) = ( ( u(k,j,i+1) - u(k,j,i) ) * rho_air(k) * ddx + & 410 431 ( v(k,j+1,i) - v(k,j,i) ) * rho_air(k) * ddy + & 411 432 ( w(k,j,i) * rho_air_zw(k) - & 412 433 w(k-1,j,i) * rho_air_zw(k-1) ) * ddzw(k) & 413 ) * ddt_3d * d_weight_pres 434 ) * ddt_3d * d_weight_pres & 435 * MERGE( 1.0_wp, 0.0_wp, & 436 BTEST( wall_flags_0(k,j,i), 0 ) & 437 ) 414 438 ENDDO 415 439 ! 416 440 !-- Compute possible PE-sum of divergences for flow_statistics 417 DO k = nzb_s_inner(j,i)+1, nzt 418 threadsum = threadsum + ABS( d(k,j,i) ) 441 DO k = nzb+1, nzt 442 threadsum = threadsum + ABS( d(k,j,i) ) & 443 * MERGE( 1.0_wp, 0.0_wp, & 444 BTEST( wall_flags_0(k,j,i), 0 ) & 445 ) 419 446 ENDDO 420 447 … … 438 465 ( w(k,j,i) * rho_air_zw(k) - & 439 466 w(k-1,j,i) * rho_air_zw(k-1) ) * ddzw(k) & 440 ) * ddt_3d * d_weight_pres * rflags_s_inner(k,j,i) 467 ) * ddt_3d * d_weight_pres & 468 * MERGE( 1.0_wp, 0.0_wp, & 469 BTEST( wall_flags_0(k,j,i), 0 ) & 470 ) 441 471 ENDDO 442 472 ENDDO … … 484 514 !-- Store computed perturbation pressure and set boundary condition in 485 515 !-- z-direction 486 !$OMP PARALLEL DO 516 !$OMP PARALLEL DO PRIVATE (i,j,k) 487 517 DO i = nxl, nxr 488 518 DO j = nys, nyn … … 499 529 IF ( ibc_p_b == 1 ) THEN 500 530 ! 501 !-- Neumann (dp/dz = 0) 502 !$OMP PARALLEL DO 503 DO i = nxlg, nxrg 504 DO j = nysg, nyng 505 tend(nzb_s_inner(j,i),j,i) = tend(nzb_s_inner(j,i)+1,j,i) 506 ENDDO 531 !-- Neumann (dp/dz = 0). Using surfae data type, first for non-natural 532 !-- surfaces, then for natural and urban surfaces 533 !-- Upward facing 534 !$OMP PARALLEL DO PRIVATE( i, j, k ) 535 DO m = 1, bc_h(0)%ns 536 i = bc_h(0)%i(m) 537 j = bc_h(0)%j(m) 538 k = bc_h(0)%k(m) 539 tend(k-1,j,i) = tend(k,j,i) 540 ENDDO 541 ! 542 !-- Downward facing 543 !$OMP PARALLEL DO PRIVATE( i, j, k ) 544 DO m = 1, bc_h(1)%ns 545 i = bc_h(1)%i(m) 546 j = bc_h(1)%j(m) 547 k = bc_h(1)%k(m) 548 tend(k+1,j,i) = tend(k,j,i) 507 549 ENDDO 508 550 509 551 ELSE 510 552 ! 511 !-- Dirichlet 512 !$OMP PARALLEL DO 513 DO i = nxlg, nxrg 514 DO j = nysg, nyng 515 tend(nzb_s_inner(j,i),j,i) = 0.0_wp 516 ENDDO 553 !-- Dirichlet. Using surface data type, first for non-natural 554 !-- surfaces, then for natural and urban surfaces 555 !-- Upward facing 556 !$OMP PARALLEL DO PRIVATE( i, j, k ) 557 DO m = 1, bc_h(0)%ns 558 i = bc_h(0)%i(m) 559 j = bc_h(0)%j(m) 560 k = bc_h(0)%k(m) 561 tend(k-1,j,i) = 0.0_wp 562 ENDDO 563 ! 564 !-- Downward facing 565 !$OMP PARALLEL DO PRIVATE( i, j, k ) 566 DO m = 1, bc_h(1)%ns 567 i = bc_h(1)%i(m) 568 j = bc_h(1)%j(m) 569 k = bc_h(1)%k(m) 570 tend(k+1,j,i) = 0.0_wp 517 571 ENDDO 518 572 … … 524 578 ! 525 579 !-- Neumann 526 !$OMP PARALLEL DO 580 !$OMP PARALLEL DO PRIVATE (i,j,k) 527 581 DO i = nxlg, nxrg 528 582 DO j = nysg, nyng … … 534 588 ! 535 589 !-- Dirichlet 536 !$OMP PARALLEL DO 590 !$OMP PARALLEL DO PRIVATE (i,j,k) 537 591 DO i = nxlg, nxrg 538 592 DO j = nysg, nyng … … 646 700 DO j = nys, nyn 647 701 648 DO k = 1, nzt 649 IF ( k > nzb_w_inner(j,i) ) THEN 650 w(k,j,i) = w(k,j,i) - dt_3d * & 651 ( tend(k+1,j,i) - tend(k,j,i) ) * ddzu(k+1) * & 652 weight_pres_l 653 ENDIF 654 ENDDO 655 656 DO k = 1, nzt 657 IF ( k > nzb_u_inner(j,i) ) THEN 658 u(k,j,i) = u(k,j,i) - dt_3d * & 659 ( tend(k,j,i) - tend(k,j,i-1) ) * ddx * & 660 weight_pres_l 661 ENDIF 662 ENDDO 663 664 DO k = 1, nzt 665 IF ( k > nzb_v_inner(j,i) ) THEN 666 v(k,j,i) = v(k,j,i) - dt_3d * & 667 ( tend(k,j,i) - tend(k,j-1,i) ) * ddy * & 668 weight_pres_l 669 ENDIF 702 DO k = nzb+1, nzt 703 w(k,j,i) = w(k,j,i) - dt_3d * & 704 ( tend(k+1,j,i) - tend(k,j,i) ) * ddzu(k+1) & 705 * weight_pres_l & 706 * MERGE( 1.0_wp, 0.0_wp, & 707 BTEST( wall_flags_0(k,j,i), 3 ) & 708 ) 709 ENDDO 710 711 DO k = nzb+1, nzt 712 u(k,j,i) = u(k,j,i) - dt_3d * & 713 ( tend(k,j,i) - tend(k,j,i-1) ) * ddx & 714 * weight_pres_l & 715 * MERGE( 1.0_wp, 0.0_wp, & 716 BTEST( wall_flags_0(k,j,i), 1 ) & 717 ) 718 ENDDO 719 720 DO k = nzb+1, nzt 721 v(k,j,i) = v(k,j,i) - dt_3d * & 722 ( tend(k,j,i) - tend(k,j-1,i) ) * ddy & 723 * weight_pres_l & 724 * MERGE( 1.0_wp, 0.0_wp, & 725 BTEST( wall_flags_0(k,j,i), 2 ) & 726 ) 670 727 ENDDO 671 728 … … 676 733 ! 677 734 !-- Sum up the volume flow through the right and north boundary 678 IF ( conserve_volume_flow .AND. bc_lr_cyc .AND. bc_ns_cyc .AND. &735 IF ( conserve_volume_flow .AND. bc_lr_cyc .AND. bc_ns_cyc .AND. & 679 736 nxr == nx ) THEN 680 737 … … 683 740 DO j = nys, nyn 684 741 !$OMP CRITICAL 685 DO k = nzb_u_inner(j,nx) + 1, nzt 686 volume_flow_l(1) = volume_flow_l(1) + u(k,j,nx) * dzw(k) 742 DO k = nzb+1, nzt 743 volume_flow_l(1) = volume_flow_l(1) + u(k,j,nxr) * dzw(k) & 744 * MERGE( 1.0_wp, 0.0_wp, & 745 BTEST( wall_flags_0(k,j,nxr), 1 )& 746 ) 687 747 ENDDO 688 748 !$OMP END CRITICAL … … 692 752 ENDIF 693 753 694 IF ( conserve_volume_flow .AND. bc_ns_cyc .AND. bc_lr_cyc .AND. &754 IF ( conserve_volume_flow .AND. bc_ns_cyc .AND. bc_lr_cyc .AND. & 695 755 nyn == ny ) THEN 696 756 … … 699 759 DO i = nxl, nxr 700 760 !$OMP CRITICAL 701 DO k = nzb_v_inner(ny,i) + 1, nzt 702 volume_flow_l(2) = volume_flow_l(2) + v(k,ny,i) * dzw(k) 761 DO k = nzb+1, nzt 762 volume_flow_l(2) = volume_flow_l(2) + v(k,nyn,i) * dzw(k) & 763 * MERGE( 1.0_wp, 0.0_wp, & 764 BTEST( wall_flags_0(k,nyn,i), 2 )& 765 ) 703 766 ENDDO 704 767 !$OMP END CRITICAL … … 727 790 DO i = nxl, nxr 728 791 DO j = nys, nyn 729 DO k = nzb_u_inner(j,i) + 1, nzt 730 u(k,j,i) = u(k,j,i) + volume_flow_offset(1) 731 ENDDO 732 DO k = nzb_v_inner(j,i) + 1, nzt 733 v(k,j,i) = v(k,j,i) + volume_flow_offset(2) 792 DO k = nzb+1, nzt 793 u(k,j,i) = u(k,j,i) + volume_flow_offset(1) & 794 * MERGE( 1.0_wp, 0.0_wp, & 795 BTEST( wall_flags_0(k,j,i), 1 ) & 796 ) 797 ENDDO 798 DO k = nzb+1, nzt 799 v(k,j,i) = v(k,j,i) + volume_flow_offset(2) & 800 * MERGE( 1.0_wp, 0.0_wp, & 801 BTEST( wall_flags_0(k,j,i), 2 ) & 802 ) 734 803 ENDDO 735 804 ENDDO … … 768 837 DO i = nxl, nxr 769 838 DO j = nys, nyn 770 DO k = nzb_s_inner(j,i)+1, nzt 771 d(k,j,i) = ( u(k,j,i+1) - u(k,j,i) ) * rho_air(k) * ddx + & 772 ( v(k,j+1,i) - v(k,j,i) ) * rho_air(k) * ddy + & 773 ( w(k,j,i) * rho_air_zw(k) - & 774 w(k-1,j,i) * rho_air_zw(k-1) ) * ddzw(k) 839 DO k = nzb+1, nzt 840 d(k,j,i) = ( ( u(k,j,i+1) - u(k,j,i) ) * rho_air(k) * ddx + & 841 ( v(k,j+1,i) - v(k,j,i) ) * rho_air(k) * ddy + & 842 ( w(k,j,i) * rho_air_zw(k) - & 843 w(k-1,j,i) * rho_air_zw(k-1) ) * ddzw(k) & 844 ) * MERGE( 1.0_wp, 0.0_wp, & 845 BTEST( wall_flags_0(k,j,i), 0 ) & 846 ) 775 847 ENDDO 776 848 DO k = nzb+1, nzt … … 783 855 DO i = nxl, nxr 784 856 DO j = nys, nyn 785 DO k = 1, nzt857 DO k = nzb+1, nzt 786 858 d(k,j,i) = ( ( u(k,j,i+1) - u(k,j,i) ) * rho_air(k) * ddx + & 787 859 ( v(k,j+1,i) - v(k,j,i) ) * rho_air(k) * ddy + & 788 860 ( w(k,j,i) * rho_air_zw(k) - & 789 861 w(k-1,j,i) * rho_air_zw(k-1) ) * ddzw(k) & 790 ) * rflags_s_inner(k,j,i) 862 ) * MERGE( 1.0_wp, 0.0_wp, & 863 BTEST( wall_flags_0(k,j,i), 0 ) & 864 ) 791 865 ENDDO 792 866 ENDDO
Note: See TracChangeset
for help on using the changeset viewer.
