Changeset 707 for palm/trunk/SOURCE/pres.f90
 Timestamp:
 Mar 29, 2011 11:39:40 AM (12 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

palm/trunk/SOURCE/pres.f90
r695 r707 4 4 ! Current revisions: 5 5 !  6 ! 6 ! Calculation of weighted average of p is now handled in the same way 7 ! regardless of the number of ghost layers (advection scheme), 8 ! multigrid and sor method are using p_loc for iterative advancements of 9 ! pressure, 10 ! localsum calculation modified for proper OpenMP reduction, 11 ! bc_lr/ns replaced by bc_lr/ns_cyc 7 12 ! 8 13 ! Former revisions: … … 99 104 100 105 ! 101 ! Multigrid method expects 1 additional grid point for the arrays102 ! d, tend and p106 ! Multigrid method expects array d to have one ghost layer. 107 ! 103 108 IF ( psolver == 'multigrid' ) THEN 104 109 105 110 DEALLOCATE( d ) 106 111 ALLOCATE( d(nzb:nzt+1,nys1:nyn+1,nxl1:nxr+1) ) 112 113 ! 114 ! Since p is later used to hold the weighted average of the substeps, it 115 ! cannot be used in the iterative solver. Therefore, its initial value is 116 ! stored on p_loc, which is then iteratively advanced in every substep. 117 IF ( intermediate_timestep_count == 1 ) THEN 118 DO i = nxl1, nxr+1 119 DO j = nys1, nyn+1 120 DO k = nzb, nzt+1 121 p_loc(k,j,i) = p(k,j,i) 122 ENDDO 123 ENDDO 124 ENDDO 125 ENDIF 107 126 108 IF ( ws_scheme_mom .OR. ws_scheme_sca ) THEN 109 110 DEALLOCATE( tend ) 111 ALLOCATE( tend(nzb:nzt+1,nys1:nyn+1,nxl1:nxr+1) ) 112 DEALLOCATE( p ) 113 ALLOCATE( p(nzb:nzt+1,nys1:nyn+1,nxl1:nxr+1) ) 114 115 ENDIF 116 127 ELSEIF ( psolver == 'sor' .AND. intermediate_timestep_count == 1 ) THEN 128 129 ! 130 ! Since p is later used to hold the weighted average of the substeps, it 131 ! cannot be used in the iterative solver. Therefore, its initial value is 132 ! stored on p_loc, which is then iteratively advanced in every substep. 133 p_loc = p 134 117 135 ENDIF 118 136 … … 297 315 ENDDO 298 316 299 localsum = ( localsum + threadsum )* dt_3d * &300 weight_pres(intermediate_timestep_count)317 localsum = localsum + threadsum * dt_3d * & 318 weight_pres(intermediate_timestep_count) 301 319 302 320 !$OMP END PARALLEL … … 362 380 ENDDO 363 381 ENDDO 364 localsum = ( localsum + threadsum ) * dt_3d&365 *weight_pres(intermediate_timestep_count)382 localsum = localsum + threadsum * dt_3d * & 383 weight_pres(intermediate_timestep_count) 366 384 !$OMP END PARALLEL 367 385 #endif … … 371 389 ! of the total domain 372 390 sums_divold_l(0:statistic_regions) = localsum 373 374 !375 ! Determine absolute minimum/maximum (only for test cases, therefore as376 ! comment line)377 ! CALL global_min_max( nzb+1, nzt, nys, nyn, nxl, nxr, d, 'abs', divmax, &378 ! divmax_ijk )379 391 380 392 CALL cpu_log( log_point_s(1), 'divergence', 'stop' ) … … 496 508 ! Solve Poisson equation for perturbation pressure using SORRed/Black 497 509 ! scheme 498 CALL sor( d, ddzu_pres, ddzw, p )499 tend = p 510 CALL sor( d, ddzu_pres, ddzw, p_loc ) 511 tend = p_loc 500 512 501 513 ELSEIF ( psolver == 'multigrid' ) THEN … … 506 518 ! to discern data exchange in multigrid ( 1 ghostpoint ) and normal grid 507 519 ! ( nbgp ghost points ). 508 exchange_mg = .TRUE.509 520 CALL poismg( tend ) 510 exchange_mg = .FALSE. 521 511 522 ! 512 523 ! Restore perturbation pressure on tend because this array is used 513 524 ! further below to correct the velocity fields 514 515 tend = p 516 IF( ws_scheme_mom .OR. ws_scheme_sca ) THEN 517 ! 518 ! Allocate p to its normal size and restore pressure. 519 DEALLOCATE( p ) 520 ALLOCATE( p(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 521 522 ENDIF 523 524 ENDIF 525 526 ! 527 ! Store perturbation pressure on array p, used in the momentum equations 528 IF ( psolver(1:7) == 'poisfft' ) THEN 529 530 IF ( intermediate_timestep_count == 1 ) THEN 531 !$OMP PARALLEL PRIVATE (i,j,k) 532 !$OMP DO 533 DO i = nxlg, nxrg 534 DO j = nysg, nyng 535 DO k = nzb, nzt+1 536 p(k,j,i) = tend(k,j,i) & 537 * weight_substep(intermediate_timestep_count) 538 ENDDO 539 ENDDO 540 ENDDO 541 !$OMP END PARALLEL 542 543 ELSE 544 !$OMP PARALLEL PRIVATE (i,j,k) 545 !$OMP DO 546 DO i = nxlg, nxrg 547 DO j = nysg, nyng 548 DO k = nzb, nzt+1 549 p(k,j,i) = p(k,j,i) + tend(k,j,i) & 550 * weight_substep(intermediate_timestep_count) 551 ENDDO 552 ENDDO 553 ENDDO 554 !$OMP END PARALLEL 555 556 ENDIF 525 DO i = nxl1, nxr+1 526 DO j = nys1, nyn+1 527 DO k = nzb, nzt+1 528 tend(k,j,i) = p_loc(k,j,i) 529 ENDDO 530 ENDDO 531 ENDDO 532 533 ENDIF 534 535 ! 536 ! Store perturbation pressure on array p, used for pressure data output. 537 ! Ghost layers are added in the output routines (except sormethod: see below) 538 IF ( intermediate_timestep_count == 1 ) THEN 539 !$OMP PARALLEL PRIVATE (i,j,k) 540 !$OMP DO 541 DO i = nxl1, nxr+1 542 DO j = nys1, nyn+1 543 DO k = nzb, nzt+1 544 p(k,j,i) = tend(k,j,i) * & 545 weight_substep(intermediate_timestep_count) 546 ENDDO 547 ENDDO 548 ENDDO 549 !$OMP END PARALLEL 550 551 ELSE 552 !$OMP PARALLEL PRIVATE (i,j,k) 553 !$OMP DO 554 DO i = nxl1, nxr+1 555 DO j = nys1, nyn+1 556 DO k = nzb, nzt+1 557 p(k,j,i) = p(k,j,i) + tend(k,j,i) * & 558 weight_substep(intermediate_timestep_count) 559 ENDDO 560 ENDDO 561 ENDDO 562 !$OMP END PARALLEL 563 564 ENDIF 557 565 558 ENDIF 566 ! 567 ! SORmethod needs ghost layers for the next timestep 568 IF ( psolver == 'sor' ) CALL exchange_horiz( p, nbgp ) 559 569 560 570 ! 561 571 ! Correction of the provisional velocities with the current perturbation 562 572 ! pressure just computed 563 IF ( conserve_volume_flow .AND. & 564 ( bc_lr == 'cyclic' .OR. bc_ns == 'cyclic' ) ) THEN 573 IF ( conserve_volume_flow .AND. ( bc_lr_cyc .OR. bc_ns_cyc ) ) THEN 565 574 volume_flow_l(1) = 0.0 566 575 volume_flow_l(2) = 0.0 567 576 ENDIF 577 568 578 !$OMP PARALLEL PRIVATE (i,j,k) 569 579 !$OMP DO … … 571 581 DO j = nys, nyn 572 582 DO k = nzb_w_inner(j,i)+1, nzt 573 w(k,j,i) = w(k,j,i)  dt_3d * 574 ( tend(k+1,j,i)  tend(k,j,i) ) * ddzu(k+1) 575 *weight_pres(intermediate_timestep_count)583 w(k,j,i) = w(k,j,i)  dt_3d * & 584 ( tend(k+1,j,i)  tend(k,j,i) ) * ddzu(k+1) * & 585 weight_pres(intermediate_timestep_count) 576 586 ENDDO 577 587 DO k = nzb_u_inner(j,i)+1, nzt 578 588 u(k,j,i) = u(k,j,i)  dt_3d * & 579 ( tend(k,j,i)  tend(k,j,i1) ) * ddx 580 *weight_pres(intermediate_timestep_count)589 ( tend(k,j,i)  tend(k,j,i1) ) * ddx * & 590 weight_pres(intermediate_timestep_count) 581 591 ENDDO 582 592 DO k = nzb_v_inner(j,i)+1, nzt 583 593 v(k,j,i) = v(k,j,i)  dt_3d * & 584 ( tend(k,j,i)  tend(k,j1,i) ) * ddy 585 *weight_pres(intermediate_timestep_count)594 ( tend(k,j,i)  tend(k,j1,i) ) * ddy * & 595 weight_pres(intermediate_timestep_count) 586 596 ENDDO 587 597 ! 588 598 ! Sum up the volume flow through the right and north boundary 589 IF ( conserve_volume_flow .AND. bc_lr == 'cyclic'.AND. &590 bc_ns == 'cyclic' .AND.i == nx ) THEN599 IF ( conserve_volume_flow .AND. bc_lr_cyc .AND. bc_ns_cyc .AND. & 600 i == nx ) THEN 591 601 !$OMP CRITICAL 592 602 DO k = nzb_2d(j,i) + 1, nzt … … 595 605 !$OMP END CRITICAL 596 606 ENDIF 597 IF ( conserve_volume_flow .AND. bc_ns == 'cyclic'.AND. &598 bc_lr == 'cyclic' .AND.j == ny ) THEN607 IF ( conserve_volume_flow .AND. bc_ns_cyc .AND. bc_lr_cyc .AND. & 608 j == ny ) THEN 599 609 !$OMP CRITICAL 600 610 DO k = nzb_2d(j,i) + 1, nzt … … 608 618 !$OMP END PARALLEL 609 619 610 IF ( psolver == 'multigrid' .OR. psolver == 'sor' ) THEN611 IF ( intermediate_timestep_count == 1 .OR. simulated_time == 0) THEN612 !$OMP PARALLEL PRIVATE (i,j,k)613 !$OMP DO614 DO i = nxl, nxr615 DO j = nys, nyn616 DO k = nzb, nzt+1617 p_sub(k,j,i) = tend(k,j,i) &618 * weight_substep(intermediate_timestep_count)619 ENDDO620 ENDDO621 ENDDO622 !$OMP END PARALLEL623 ELSE624 !$OMP PARALLEL PRIVATE (i,j,k)625 !$OMP DO626 DO i = nxl, nxr627 DO j = nys, nyn628 DO k = nzb, nzt+1629 p_sub(k,j,i) = p_sub(k,j,i) + tend(k,j,i) &630 * weight_substep(intermediate_timestep_count)631 ENDDO632 ENDDO633 ENDDO634 !$OMP END PARALLEL635 ENDIF636 637 IF ( intermediate_timestep_count == intermediate_timestep_count_max ) &638 THEN639 !$OMP PARALLEL PRIVATE (i,j,k)640 !$OMP DO641 DO i = nxl, nxr642 DO j = nys, nyn643 DO k = nzb, nzt+1644 p(k,j,i) = p_sub(k,j,i)645 ENDDO646 ENDDO647 ENDDO648 !$OMP END PARALLEL649 ENDIF650 ENDIF651 652 !653 ! Resize tend to its normal size in case of multigrid and wsscheme.654 IF ( psolver == 'multigrid' .AND. ( ws_scheme_mom &655 .OR. ws_scheme_sca ) ) THEN656 DEALLOCATE( tend )657 ALLOCATE( tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )658 ENDIF659 660 661 620 ! 662 621 ! Conserve the volume flow 663 IF ( conserve_volume_flow .AND. & 664 ( bc_lr == 'cyclic' .AND. bc_ns == 'cyclic' ) ) THEN 622 IF ( conserve_volume_flow .AND. ( bc_lr_cyc .AND. bc_ns_cyc ) ) THEN 665 623 666 624 #if defined( __parallel )
Note: See TracChangeset
for help on using the changeset viewer.