- Timestamp:
- Mar 13, 2014 2:30:59 PM (11 years ago)
- Location:
- palm/trunk/SOURCE
- Files:
-
- 1 added
- 3 edited
-
Makefile (modified) (2 diffs)
-
mod_kinds.f90 (added)
-
poisfft.f90 (modified) (9 diffs)
-
pres.f90 (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/Makefile
r1242 r1306 20 20 # Current revisions: 21 21 # ------------------ 22 # 22 # +mod_kinds 23 23 # 24 24 # Former revisions: … … 161 161 lpm_set_attributes.f90 lpm_sort_arrays.f90 \ 162 162 lpm_write_exchange_statistics.f90 lpm_write_restart_file.f90 \ 163 ls_forcing.f90 \164 message.f90 microphysics.f90 modules.f90netcdf.f90 nudging.f90 \165 163 ls_forcing.f90 message.f90 microphysics.f90 modules.f90 mod_kinds.f90 \ 164 netcdf.f90 nudging.f90 \ 165 package_parin.f90 palm.f90 parin.f90 plant_canopy_model.f90 poisfft.f90 \ 166 166 poismg.f90 prandtl_fluxes.f90 pres.f90 print_1d.f90 \ 167 167 production_e.f90 prognostic_equations.f90 random_function.f90 \ -
palm/trunk/SOURCE/poisfft.f90
r1217 r1306 15 15 ! PALM. If not, see <http://www.gnu.org/licenses/>. 16 16 ! 17 ! Copyright 1997-201 2Leibniz University Hannover17 ! Copyright 1997-2014 Leibniz University Hannover 18 18 !--------------------------------------------------------------------------------! 19 19 ! 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! 22 ! openmp sections removed from the overlap branch, 23 ! second argument removed from parameter list 23 24 ! 24 25 ! Former revisions: … … 146 147 ! Description: 147 148 ! ------------ 148 ! See below. 149 ! Original version by Stephan Siano (pois3d), as of July 23, 1996 150 ! Adapted for 2D-domain-decomposition by Siegfried Raasch, July 3, 1997 151 ! 152 ! Solves the Poisson equation with a 2D spectral method 153 ! d^2 p / dx^2 + d^2 p / dy^2 + d^2 p / dz^2 = s 154 ! 155 ! Input: 156 ! real ar contains (nnz,nny,nnx) elements of the velocity divergence, 157 ! starting from (1,nys,nxl) 158 ! 159 ! Output: 160 ! real ar contains the solution for perturbation pressure p 149 161 !------------------------------------------------------------------------------! 150 151 !--------------------------------------------------------------------------!152 ! poisfft !153 ! !154 ! Original version: Stephan Siano (pois3d) !155 ! !156 ! Institute of Meteorology and Climatology, University of Hannover !157 ! Germany !158 ! !159 ! Version as of July 23,1996 !160 ! !161 ! !162 ! Version for parallel computers: Siegfried Raasch !163 ! !164 ! Version as of July 03,1997 !165 ! !166 ! Solves the Poisson equation with a 2D spectral method !167 ! d^2 p / dx^2 + d^2 p / dy^2 + d^2 p / dz^2 = s !168 ! !169 ! Input: !170 ! real ar contains in the (nnx,nny,nnz) elements, !171 ! starting from the element (1,nys,nxl), the !172 ! values for s !173 ! real work Temporary array !174 ! !175 ! Output: !176 ! real ar contains the solution for p !177 !--------------------------------------------------------------------------!178 162 179 163 USE fft_xy … … 227 211 228 212 #if ! defined ( __check ) 229 SUBROUTINE poisfft( ar , ar_inv_test)213 SUBROUTINE poisfft( ar ) 230 214 231 215 USE control_parameters, ONLY : fft_method, transpose_compute_overlap … … 236 220 IMPLICIT NONE 237 221 238 INTEGER :: ind_even, ind_odd, ind_third, ii, iind, inew, jj, jind, & 239 jnew, ki, kk, knew, n, nblk, nnx_y, nny_z, nnz_t, nnz_x, & 240 nxl_y_bound, nxr_y_bound 222 INTEGER :: ii, iind, inew, jj, jind, jnew, ki, kk, knew, n, nblk, & 223 nnx_y, nny_z, nnz_t, nnz_x, nxl_y_bound, nxr_y_bound 241 224 INTEGER, DIMENSION(4) :: isave 242 225 243 226 REAL, DIMENSION(1:nz,nys:nyn,nxl:nxr) :: ar 244 REAL, DIMENSION(nys:nyn,nxl:nxr,1:nz) :: ar_inv_test ! work array tend from pres245 227 !$acc declare create( ar_inv ) 246 228 REAL, DIMENSION(nys:nyn,nxl:nxr,1:nz) :: ar_inv 247 229 248 REAL, DIMENSION(:,:,:), ALLOCATABLE :: f_in, f_inv, f_out_y, f_out_z249 REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: ar1230 REAL, DIMENSION(:,:,:), ALLOCATABLE :: ar1, f_in, f_inv, f_out_y, & 231 f_out_z 250 232 251 233 … … 408 390 sendrecvcount_zx = nnx * nny * nnz_x 409 391 410 ALLOCATE( ar1(0:nx,nys_x:nyn_x,nzb_x:nzt_x ,2) )392 ALLOCATE( ar1(0:nx,nys_x:nyn_x,nzb_x:nzt_x) ) 411 393 ALLOCATE( f_in(nys:nyn,nxl:nxr,1:nz) ) 412 394 413 DO kk = 1, nblk+1 414 ind_odd = MOD( kk, 2 ) + 1 415 ind_even = MOD( kk+1, 2 ) + 1 416 !$OMP sections private(ki,knew,n) 417 !$OMP section 418 IF ( kk <= nblk ) THEN 419 420 IF ( kk == 1 ) THEN 421 CALL cpu_log( log_point_s(5), 'transpo forward', 'start' ) 422 ELSE 423 CALL cpu_log( log_point_s(5), 'transpo forward', 'continue' ) 424 ENDIF 425 426 DO knew = 1, nz 427 ki = kk + nblk * ( knew - 1 ) 428 f_in(:,:,knew) = f_inv(:,:,ki) 429 ENDDO 430 431 CALL transpose_zx( f_in, ar1(:,:,:,ind_odd)) 432 CALL cpu_log( log_point_s(5), 'transpo forward', 'pause' ) 433 395 DO kk = 1, nblk 396 397 IF ( kk == 1 ) THEN 398 CALL cpu_log( log_point_s(5), 'transpo forward', 'start' ) 399 ELSE 400 CALL cpu_log( log_point_s(5), 'transpo forward', 'continue' ) 434 401 ENDIF 435 402 436 !$OMP section 437 IF ( kk >= 2 ) THEN 438 439 IF ( kk == 2 ) THEN 440 CALL cpu_log( log_point_s(4), 'fft_x', 'start' ) 441 ELSE 442 CALL cpu_log( log_point_s(4), 'fft_x', 'continue' ) 443 ENDIF 444 445 n = isave(2) + kk - 2 446 CALL fft_x( ar1(:,:,:,ind_even), 'forward', ar_2d = f_out_z(:,:,n)) 447 CALL cpu_log( log_point_s(4), 'fft_x', 'pause' ) 448 403 DO knew = 1, nz 404 ki = kk + nblk * ( knew - 1 ) 405 f_in(:,:,knew) = f_inv(:,:,ki) 406 ENDDO 407 408 CALL transpose_zx( f_in, ar1(:,:,:)) 409 CALL cpu_log( log_point_s(5), 'transpo forward', 'pause' ) 410 411 IF ( kk == 1 ) THEN 412 CALL cpu_log( log_point_s(4), 'fft_x', 'start' ) 413 ELSE 414 CALL cpu_log( log_point_s(4), 'fft_x', 'continue' ) 449 415 ENDIF 450 !$OMP end sections 416 417 n = isave(2) + kk - 1 418 CALL fft_x( ar1(:,:,:), 'forward', ar_2d = f_out_z(:,:,n)) 419 CALL cpu_log( log_point_s(4), 'fft_x', 'pause' ) 451 420 452 421 ENDDO … … 479 448 sendrecvcount_xy = nnx_y * ( nyn_x-nys_x+1 ) * ( nzt_x-nzb_x+1 ) 480 449 481 ALLOCATE( ar1(0:ny,nxl_y:nxr_y,nzb_y:nzt_y ,2) )450 ALLOCATE( ar1(0:ny,nxl_y:nxr_y,nzb_y:nzt_y) ) 482 451 ALLOCATE( f_in(nys_x:nyn_x,nzb_x:nzt_x,0:nx) ) 483 452 484 DO ii = 0, nblk+1 485 ind_odd = MOD( ii+1, 2 ) + 1 486 ind_even = MOD( ii+2, 2 ) + 1 487 !$OMP sections private(ki,knew,n) 488 !$OMP section 489 IF ( ii <= nblk ) THEN 490 491 CALL cpu_log( log_point_s(5), 'transpo forward', 'continue' ) 492 493 DO inew = 0, nx-1 494 iind = ii + ( nblk + 1 ) * inew 495 f_in(:,:,inew) = f_inv(:,:,iind) 496 ENDDO 497 498 CALL transpose_xy( f_in, ar1(:,:,:,ind_odd) ) 499 500 CALL cpu_log( log_point_s(5), 'transpo forward', 'pause' ) 501 453 DO ii = 0, nblk 454 455 CALL cpu_log( log_point_s(5), 'transpo forward', 'continue' ) 456 457 DO inew = 0, nx-1 458 iind = ii + ( nblk + 1 ) * inew 459 f_in(:,:,inew) = f_inv(:,:,iind) 460 ENDDO 461 462 CALL transpose_xy( f_in, ar1(:,:,:) ) 463 464 CALL cpu_log( log_point_s(5), 'transpo forward', 'pause' ) 465 466 IF ( ii == 1 ) THEN 467 CALL cpu_log( log_point_s(7), 'fft_y', 'start' ) 468 ELSE 469 CALL cpu_log( log_point_s(7), 'fft_y', 'continue' ) 502 470 ENDIF 503 471 504 !$OMP section 505 IF ( ii >= 1 ) THEN 506 507 IF ( ii == 1 ) THEN 508 CALL cpu_log( log_point_s(7), 'fft_y', 'start' ) 509 ELSE 510 CALL cpu_log( log_point_s(7), 'fft_y', 'continue' ) 511 ENDIF 512 513 nxl_y_bound = isave(2) 514 nxr_y_bound = isave(3) 515 n = isave(2) + ii - 1 516 ! CALL fft_y( ar1(:,:,:,ind_even), 'forward', ar_3d = f_out_y, & 517 ! ni = n ) 518 CALL fft_y( ar1(:,:,:,ind_even), 'forward', ar_tr = f_out_y, & 519 nxl_y_bound = nxl_y_bound, nxr_y_bound = nxr_y_bound, & 520 nxl_y_l = n, nxr_y_l = n ) 521 522 CALL cpu_log( log_point_s(7), 'fft_y', 'pause' ) 523 524 ENDIF 525 !$OMP end sections 472 nxl_y_bound = isave(2) 473 nxr_y_bound = isave(3) 474 n = isave(2) + ii 475 CALL fft_y( ar1(:,:,:), 'forward', ar_tr = f_out_y, & 476 nxl_y_bound = nxl_y_bound, nxr_y_bound = nxr_y_bound, & 477 nxl_y_l = n, nxr_y_l = n ) 478 479 CALL cpu_log( log_point_s(7), 'fft_y', 'pause' ) 526 480 527 481 ENDDO … … 554 508 sendrecvcount_yz = ( nxr_y-nxl_y+1 ) * nny_z * ( nzt_y-nzb_y+1 ) 555 509 556 ALLOCATE( ar1(nxl_z:nxr_z,nys_z:nyn_z,1:nz ,3) )510 ALLOCATE( ar1(nxl_z:nxr_z,nys_z:nyn_z,1:nz) ) 557 511 ALLOCATE( f_in(nxl_y:nxr_y,nzb_y:nzt_y,0:ny) ) 558 512 559 DO jj = 0, nblk+2 560 ind_odd = MOD( jj+3, 3 ) + 1 561 ind_even = MOD( jj+2, 3 ) + 1 562 ind_third = MOD( jj+1, 3 ) + 1 563 !$OMP sections private(ki,knew,n) 564 !$OMP section 565 IF ( jj <= nblk ) THEN 566 ! 567 !-- Forward Fourier Transformation 568 !-- Transposition y --> z 569 CALL cpu_log( log_point_s(5), 'transpo forward', 'continue' ) 570 571 DO jnew = 0, ny-1 572 jind = jj + ( nblk + 1 ) * jnew 573 f_in(:,:,jnew) =f_inv(:,:,jind) 574 ENDDO 575 576 CALL transpose_yz( f_in, ar1(:,:,:,ind_odd) ) 577 578 IF ( jj == nblk ) THEN 579 CALL cpu_log( log_point_s(5), 'transpo forward', 'stop' ) 580 ELSE 581 CALL cpu_log( log_point_s(5), 'transpo forward', 'pause' ) 582 ENDIF 583 513 DO jj = 0, nblk 514 ! 515 !-- Forward Fourier Transformation 516 !-- Transposition y --> z 517 CALL cpu_log( log_point_s(5), 'transpo forward', 'continue' ) 518 519 DO jnew = 0, ny-1 520 jind = jj + ( nblk + 1 ) * jnew 521 f_in(:,:,jnew) = f_inv(:,:,jind) 522 ENDDO 523 524 CALL transpose_yz( f_in, ar1(:,:,:) ) 525 526 IF ( jj == nblk ) THEN 527 CALL cpu_log( log_point_s(5), 'transpo forward', 'stop' ) 528 ELSE 529 CALL cpu_log( log_point_s(5), 'transpo forward', 'pause' ) 584 530 ENDIF 585 531 586 IF ( jj >= 2 ) THEN 587 ! 588 !-- Inverse Fourier Transformation 589 !-- Transposition z --> y 590 !-- Only one thread should call MPI routines, therefore forward and 591 !-- backward tranpose are in the same section 592 IF ( jj == 2 ) THEN 593 CALL cpu_log( log_point_s(8), 'transpo invers', 'start' ) 594 ELSE 595 CALL cpu_log( log_point_s(8), 'transpo invers', 'continue' ) 596 ENDIF 597 598 CALL transpose_zy( ar1(:,:,:,ind_third), f_in ) 599 600 DO jnew = 0, ny-1 601 jind = jj-2 + ( nblk + 1 ) * jnew 602 f_inv(:,:,jind) = f_in(:,:,jnew) 603 ENDDO 604 605 CALL cpu_log( log_point_s(8), 'transpo invers', 'pause' ) 606 532 ! 533 !-- Solve the tridiagonal equation system along z 534 CALL cpu_log( log_point_s(6), 'tridia', 'start' ) 535 536 n = isave(2) + jj 537 CALL tridia_substi_overlap( ar1(:,:,:), n ) 538 539 CALL cpu_log( log_point_s(6), 'tridia', 'stop' ) 540 541 ! 542 !-- Inverse Fourier Transformation 543 !-- Transposition z --> y 544 !-- Only one thread should call MPI routines, therefore forward and 545 !-- backward tranpose are in the same section 546 IF ( jj == 0 ) THEN 547 CALL cpu_log( log_point_s(8), 'transpo invers', 'start' ) 548 ELSE 549 CALL cpu_log( log_point_s(8), 'transpo invers', 'continue' ) 607 550 ENDIF 608 551 609 !$OMP section 610 IF ( jj >= 1 .AND. jj <= nblk+1 ) THEN 611 ! 612 !-- Solve the tridiagonal equation system along z 613 CALL cpu_log( log_point_s(6), 'tridia', 'start' ) 614 615 n = isave(2) + jj - 1 616 CALL tridia_substi_overlap( ar1(:,:,:,ind_even), n ) 617 618 CALL cpu_log( log_point_s(6), 'tridia', 'stop' ) 619 ENDIF 620 !$OMP end sections 552 CALL transpose_zy( ar1(:,:,:), f_in ) 553 554 DO jnew = 0, ny-1 555 jind = jj + ( nblk + 1 ) * jnew 556 f_inv(:,:,jind) = f_in(:,:,jnew) 557 ENDDO 558 559 CALL cpu_log( log_point_s(8), 'transpo invers', 'pause' ) 621 560 622 561 ENDDO … … 650 589 sendrecvcount_xy = nnx_y * ( nyn_x-nys_x+1 ) * ( nzt_x-nzb_x+1 ) 651 590 652 ALLOCATE( ar1(0:ny,nxl_y:nxr_y,nzb_y:nzt_y ,2) )591 ALLOCATE( ar1(0:ny,nxl_y:nxr_y,nzb_y:nzt_y) ) 653 592 ALLOCATE( f_in(nys_x:nyn_x,nzb_x:nzt_x,0:nx) ) 654 593 655 DO ii = 0, nblk+1 656 ind_odd = MOD( ii+1, 2 ) + 1 657 ind_even = MOD( ii+2, 2 ) + 1 658 !$OMP sections private(ki,knew,n) 659 !$OMP section 660 IF ( ii <= nblk ) THEN 661 662 CALL cpu_log( log_point_s(7), 'fft_y', 'continue' ) 663 664 n = isave(2) + ii 665 nxl_y_bound = isave(2) 666 nxr_y_bound = isave(3) 667 668 ! CALL fft_y( ar1(:,:,:,ind_even), 'backward', ar_3d = f_out_y, & 669 ! ni = n ) 670 CALL fft_y( ar1(:,:,:,ind_even), 'backward', ar_tr = f_out_y, & 671 nxl_y_bound = nxl_y_bound, nxr_y_bound = nxr_y_bound, & 672 nxl_y_l = n, nxr_y_l = n ) 673 674 IF ( ii == nblk ) THEN 675 CALL cpu_log( log_point_s(7), 'fft_y', 'stop' ) 676 ELSE 677 CALL cpu_log( log_point_s(7), 'fft_y', 'pause' ) 678 ENDIF 679 594 DO ii = 0, nblk 595 596 CALL cpu_log( log_point_s(7), 'fft_y', 'continue' ) 597 598 n = isave(2) + ii 599 nxl_y_bound = isave(2) 600 nxr_y_bound = isave(3) 601 602 CALL fft_y( ar1(:,:,:), 'backward', ar_tr = f_out_y, & 603 nxl_y_bound = nxl_y_bound, nxr_y_bound = nxr_y_bound, & 604 nxl_y_l = n, nxr_y_l = n ) 605 606 IF ( ii == nblk ) THEN 607 CALL cpu_log( log_point_s(7), 'fft_y', 'stop' ) 608 ELSE 609 CALL cpu_log( log_point_s(7), 'fft_y', 'pause' ) 680 610 ENDIF 681 611 682 !$OMP section 683 IF ( ii >= 1 ) THEN 684 685 CALL cpu_log( log_point_s(8), 'transpo invers', 'continue' ) 686 687 CALL transpose_yx( ar1(:,:,:,ind_odd), f_in ) 688 689 DO inew = 0, nx-1 690 iind = ii-1 + (nblk+1) * inew 691 f_inv(:,:,iind) = f_in(:,:,inew) 692 ENDDO 693 694 CALL cpu_log( log_point_s(8), 'transpo invers', 'pause' ) 695 696 ENDIF 697 !$OMP end sections 612 CALL cpu_log( log_point_s(8), 'transpo invers', 'continue' ) 613 614 CALL transpose_yx( ar1(:,:,:), f_in ) 615 616 DO inew = 0, nx-1 617 iind = ii + (nblk+1) * inew 618 f_inv(:,:,iind) = f_in(:,:,inew) 619 ENDDO 620 621 CALL cpu_log( log_point_s(8), 'transpo invers', 'pause' ) 698 622 699 623 ENDDO … … 727 651 sendrecvcount_zx = nnx * nny * nnz_x 728 652 729 ALLOCATE( ar1(0:nx,nys_x:nyn_x,nzb_x:nzt_x ,2) )653 ALLOCATE( ar1(0:nx,nys_x:nyn_x,nzb_x:nzt_x) ) 730 654 ALLOCATE( f_in(nys:nyn,nxl:nxr,1:nz) ) 731 655 732 DO kk = 1, nblk+1 733 ind_odd = MOD( kk, 2 ) + 1 734 ind_even = MOD( kk+1, 2 ) + 1 735 !$OMP sections private(ki,knew,n) 736 !$OMP section 737 IF ( kk <= nblk ) THEN 738 739 CALL cpu_log( log_point_s(4), 'fft_x', 'continue' ) 740 741 n = isave(2) + kk - 1 742 CALL fft_x( ar1(:,:,:,ind_even), 'backward', f_out_z(:,:,n)) 743 744 IF ( kk == nblk ) THEN 745 CALL cpu_log( log_point_s(4), 'fft_x', 'stop' ) 746 ELSE 747 CALL cpu_log( log_point_s(4), 'fft_x', 'pause' ) 748 ENDIF 749 656 DO kk = 1, nblk 657 658 CALL cpu_log( log_point_s(4), 'fft_x', 'continue' ) 659 660 n = isave(2) + kk - 1 661 CALL fft_x( ar1(:,:,:), 'backward', f_out_z(:,:,n)) 662 663 IF ( kk == nblk ) THEN 664 CALL cpu_log( log_point_s(4), 'fft_x', 'stop' ) 665 ELSE 666 CALL cpu_log( log_point_s(4), 'fft_x', 'pause' ) 750 667 ENDIF 751 668 752 !$OMP section 753 IF ( kk >= 2 ) THEN 754 755 CALL cpu_log( log_point_s(8), 'transpo invers', 'continue' ) 756 757 CALL transpose_xz( ar1(:,:,:,ind_odd), f_in ) 758 759 DO knew = 1, nz 760 ki = kk-1 + nblk * (knew-1) 761 f_inv(:,:,ki) = f_in(:,:,knew) 762 ENDDO 763 764 IF ( kk == nblk+1 ) THEN 765 CALL cpu_log( log_point_s(8), 'transpo invers', 'stop' ) 766 ELSE 767 CALL cpu_log( log_point_s(8), 'transpo invers', 'pause' ) 768 ENDIF 769 669 CALL cpu_log( log_point_s(8), 'transpo invers', 'continue' ) 670 671 CALL transpose_xz( ar1(:,:,:), f_in ) 672 673 DO knew = 1, nz 674 ki = kk + nblk * (knew-1) 675 f_inv(:,:,ki) = f_in(:,:,knew) 676 ENDDO 677 678 IF ( kk == nblk ) THEN 679 CALL cpu_log( log_point_s(8), 'transpo invers', 'stop' ) 680 ELSE 681 CALL cpu_log( log_point_s(8), 'transpo invers', 'pause' ) 770 682 ENDIF 771 !$OMP end sections772 683 773 684 ENDDO -
palm/trunk/SOURCE/pres.f90
r1258 r1306 20 20 ! Current revisions: 21 21 ! ------------------ 22 ! 22 ! second argument removed from routine poisfft 23 23 ! 24 24 ! Former revisions: … … 414 414 IF ( psolver == 'poisfft' ) THEN 415 415 416 CALL poisfft( d , tend)416 CALL poisfft( d ) 417 417 418 418 ENDIF
Note: See TracChangeset
for help on using the changeset viewer.
