Changeset 667 for palm/trunk/SOURCE/pres.f90
 Timestamp:
 Dec 23, 2010 12:06:00 PM (11 years ago)
 Location:
 palm/trunk/SOURCE
 Files:

 2 edited
Legend:
 Unmodified
 Added
 Removed

palm/trunk/SOURCE

Property
svn:mergeinfo
set to
(toggle deleted branches)
/palm/branches/suehring 423666 /palm/branches/letzel/masked_output/SOURCE 296409

Property
svn:mergeinfo
set to
(toggle deleted branches)

palm/trunk/SOURCE/pres.f90
r623 r667 4 4 ! Current revisions: 5 5 !  6 ! 6 ! New allocation of tend when wsscheme and multigrid is used. This is due to 7 ! reasons of perforance of the data_exchange. The same is done with p after 8 ! poismg is called. 9 ! nxl1, nxr+1, nys1, nyn+1 replaced by nxlg, nxrg, nysg, nyng when no 10 ! multigrid is used. Calls of exchange_horiz are modified. 11 ! 12 ! bugfix: After pressure correction no volume flow correction in case of 13 ! noncyclic boundary conditions 14 ! (has to be done only before pressure correction) 15 ! 16 ! Call of SOR routine is referenced with ddzu_pres. 7 17 ! 8 18 ! Former revisions: … … 74 84 75 85 ! 76 ! Multigrid method needs additional grid points for the divergence array 86 ! Multigrid method expects 1 additional grid point for the arrays 87 ! d, tend and p 77 88 IF ( psolver == 'multigrid' ) THEN 89 78 90 DEALLOCATE( d ) 79 ALLOCATE( d(nzb:nzt+1,nys1:nyn+1,nxl1:nxr+1) ) 91 ALLOCATE( d(nzb:nzt+1,nys1:nyn+1,nxl1:nxr+1) ) 92 93 IF ( ws_scheme_mom .OR. ws_scheme_sca ) THEN 94 95 DEALLOCATE( tend ) 96 ALLOCATE( tend(nzb:nzt+1,nys1:nyn+1,nxl1:nxr+1) ) 97 DEALLOCATE( p ) 98 ALLOCATE( p(nzb:nzt+1,nys1:nyn+1,nxl1:nxr+1) ) 99 100 ENDIF 101 80 102 ENDIF 81 103 … … 103 125 ! Sum up the volume flow through the south/north boundary 104 126 DO k = nzb_2d(j,i) + 1, nzt 105 volume_flow_l(1) = volume_flow_l(1) + u(k,j,i) * dz u(k)127 volume_flow_l(1) = volume_flow_l(1) + u(k,j,i) * dzw(k) 106 128 ENDDO 107 129 ENDDO … … 117 139 / volume_flow_area(1) 118 140 119 DO j = nys 1, nyn+1120 DO k = nzb_ v_inner(j,i) + 1, nzt141 DO j = nysg, nyng 142 DO k = nzb_2d(j,i) + 1, nzt 121 143 u(k,j,i) = u(k,j,i) + volume_flow_offset(1) 122 144 ENDDO … … 142 164 ! Sum up the volume flow through the south/north boundary 143 165 DO k = nzb_2d(j,i) + 1, nzt 144 volume_flow_l(2) = volume_flow_l(2) + v(k,j,i) * dz u(k)166 volume_flow_l(2) = volume_flow_l(2) + v(k,j,i) * dzw(k) 145 167 ENDDO 146 168 ENDDO … … 156 178 / volume_flow_area(2) 157 179 158 DO i = nxl 1, nxr+1180 DO i = nxlg, nxrg 159 181 DO k = nzb_v_inner(j,i) + 1, nzt 160 182 v(k,j,i) = v(k,j,i) + volume_flow_offset(2) … … 186 208 w_l(k) = w_l(k) / ngp_2dh_outer(k,0) 187 209 ENDDO 188 DO i = nxl 1, nxr+1189 DO j = nys 1, nyn+1210 DO i = nxlg, nxrg 211 DO j = nysg, nyng 190 212 DO k = nzb_w_inner(j,i)+1, nzt 191 213 w(k,j,i) = w(k,j,i)  w_l(k) … … 267 289 DO j = nys, nyn 268 290 DO k = nzb_s_inner(j,i)+1, nzt 269 270 271 291 d(k,j,i) = ( ( u(k,j,i+1)  u(k,j,i) ) * ddx + & 292 ( v(k,j+1,i)  v(k,j,i) ) * ddy + & 293 ( w(k,j,i)  w(k1,j,i) ) * ddzw(k) ) * ddt_3d 272 294 ENDDO 273 295 ENDDO … … 298 320 DO k = nzb_s_inner(j,i)+1, nzt 299 321 d(k,j,i) = ( ( u(k,j,i+1)  u(k,j,i) ) * ddx + & 300 301 322 ( v(k,j+1,i)  v(k,j,i) ) * ddy + & 323 ( w(k,j,i)  w(k1,j,i) ) * ddzw(k) ) * ddt_3d 302 324 ENDDO 303 325 ENDDO … … 331 353 ! comment line) 332 354 ! CALL global_min_max( nzb+1, nzt, nys, nyn, nxl, nxr, d, 'abs', divmax, & 333 ! 355 ! divmax_ijk ) 334 356 335 357 CALL cpu_log( log_point_s(1), 'divergence', 'stop' ) … … 364 386 IF ( nxra > nxr .OR. nyna > nyn .OR. nza > nz ) THEN 365 387 DEALLOCATE( tend ) 366 ALLOCATE( tend(nzb:nzt+1,nys 1:nyn+1,nxl1:nxr+1) )388 ALLOCATE( tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 367 389 ENDIF 368 390 … … 387 409 ! Neumann (dp/dz = 0) 388 410 !$OMP PARALLEL DO 389 DO i = nxl 1, nxr+1390 DO j = nys 1, nyn+1411 DO i = nxlg, nxrg 412 DO j = nysg, nyng 391 413 tend(nzb_s_inner(j,i),j,i) = tend(nzb_s_inner(j,i)+1,j,i) 392 414 ENDDO … … 400 422 ! the computation (cf. above: computation of divergences). 401 423 !$OMP PARALLEL DO 402 DO i = nxl 1, nxr+1403 DO j = nys 1, nyn+1424 DO i = nxlg, nxrg 425 DO j = nysg, nyng 404 426 tend(nzb_s_inner(j,i),j,i) = tend(nzb_s_inner(j,i)+1,j,i) 405 427 ENDDO … … 410 432 ! Dirichlet 411 433 !$OMP PARALLEL DO 412 DO i = nxl 1, nxr+1413 DO j = nys 1, nyn+1434 DO i = nxlg, nxrg 435 DO j = nysg, nyng 414 436 tend(nzb_s_inner(j,i),j,i) = 0.0 415 437 ENDDO … … 424 446 ! Neumann 425 447 !$OMP PARALLEL DO 426 DO i = nxl 1, nxr+1427 DO j = nys 1, nyn+1448 DO i = nxlg, nxrg 449 DO j = nysg, nyng 428 450 tend(nzt+1,j,i) = tend(nzt,j,i) 429 451 ENDDO … … 434 456 ! Dirichlet 435 457 !$OMP PARALLEL DO 436 DO i = nxl 1, nxr+1437 DO j = nys 1, nyn+1458 DO i = nxlg, nxrg 459 DO j = nysg, nyng 438 460 tend(nzt+1,j,i) = 0.0 439 461 ENDDO … … 444 466 ! 445 467 ! Exchange boundaries for p 446 CALL exchange_horiz( tend )468 CALL exchange_horiz( tend, nbgp ) 447 469 448 470 ELSEIF ( psolver == 'sor' ) THEN … … 451 473 ! Solve Poisson equation for perturbation pressure using SORRed/Black 452 474 ! scheme 453 CALL sor( d, ddzu , ddzw, p )475 CALL sor( d, ddzu_pres, ddzw, p ) 454 476 tend = p 455 477 … … 458 480 ! 459 481 ! Solve Poisson equation for perturbation pressure using Multigrid scheme, 460 ! array tend is used to store the residuals 482 ! array tend is used to store the residuals, logical exchange_mg is used 483 ! to discern data exchange in multigrid ( 1 ghostpoint ) and normal grid 484 ! ( nbgp ghost points ). 485 exchange_mg = .TRUE. 461 486 CALL poismg( tend ) 462 487 exchange_mg = .FALSE. 463 488 ! 464 489 ! Restore perturbation pressure on tend because this array is used 465 490 ! further below to correct the velocity fields 491 466 492 tend = p 493 IF( ws_scheme_mom .OR. ws_scheme_sca ) THEN 494 ! 495 ! Allocate p to its normal size and restore pressure. 496 DEALLOCATE( p ) 497 ALLOCATE( p(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 498 DO i = nxl, nxr 499 DO j = nys, nyn 500 DO k = nzb_s_inner(j,i), nzt 501 p(k,j,i) = tend(k,j,i) 502 ENDDO 503 ENDDO 504 ENDDO 505 ENDIF 467 506 468 507 ENDIF … … 476 515 ! optimization 477 516 !$OMP PARALLEL DO 478 DO j = nys 1, nyn+1517 DO j = nysg, nyng 479 518 DO k = nzb, nzt+1 480 p(k,j,nxl 1) = tend(k,j,nxl1)481 p(k,j,nxr+1 ) = tend(k,j,nxr+1)519 p(k,j,nxlg:nxl1) = tend(k,j,nxlg:nxl1) 520 p(k,j,nxr+1:nxrg) = tend(k,j,nxr+1:nxrg) 482 521 ENDDO 483 522 ENDDO … … 496 535 DO i = nxl, nxr 497 536 IF ( psolver(1:7) == 'poisfft' ) THEN 498 DO j = nys 1, nyn+1537 DO j = nysg, nyng 499 538 DO k = nzb, nzt+1 500 539 p(k,j,i) = tend(k,j,i) … … 517 556 ! Sum up the volume flow through the right and north boundary 518 557 IF ( conserve_volume_flow .AND. bc_lr == 'cyclic' .AND. & 519 i == nx ) THEN558 bc_ns == 'cyclic' .AND. i == nx ) THEN 520 559 !$OMP CRITICAL 521 560 DO k = nzb_2d(j,i) + 1, nzt 522 volume_flow_l(1) = volume_flow_l(1) + u(k,j,i) * dz u(k)561 volume_flow_l(1) = volume_flow_l(1) + u(k,j,i) * dzw(k) 523 562 ENDDO 524 563 !$OMP END CRITICAL 525 564 ENDIF 526 565 IF ( conserve_volume_flow .AND. bc_ns == 'cyclic' .AND. & 527 j == ny ) THEN566 bc_lr == 'cyclic' .AND. j == ny ) THEN 528 567 !$OMP CRITICAL 529 568 DO k = nzb_2d(j,i) + 1, nzt 530 volume_flow_l(2) = volume_flow_l(2) + v(k,j,i) * dz u(k)569 volume_flow_l(2) = volume_flow_l(2) + v(k,j,i) * dzw(k) 531 570 ENDDO 532 571 !$OMP END CRITICAL … … 538 577 539 578 ! 579 ! Resize tend to its normal size in case of multigrid and wsscheme. 580 IF ( psolver == 'multigrid' .AND. ( ws_scheme_mom & 581 .OR. ws_scheme_sca ) ) THEN 582 DEALLOCATE( tend ) 583 ALLOCATE( tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 584 ENDIF 585 586 ! 540 587 ! Conserve the volume flow 541 588 IF ( conserve_volume_flow .AND. & 542 ( bc_lr == 'cyclic' . OR. bc_ns == 'cyclic' ) ) THEN589 ( bc_lr == 'cyclic' .AND. bc_ns == 'cyclic' ) ) THEN 543 590 544 591 #if defined( __parallel ) … … 557 604 DO i = nxl, nxr 558 605 DO j = nys, nyn 559 IF ( bc_lr == 'cyclic' ) THEN 560 DO k = nzb_u_inner(j,i) + 1, nzt 561 u(k,j,i) = u(k,j,i) + volume_flow_offset(1) 562 ENDDO 563 ENDIF 564 IF ( bc_ns == 'cyclic' ) THEN 565 DO k = nzb_v_inner(j,i) + 1, nzt 566 v(k,j,i) = v(k,j,i) + volume_flow_offset(2) 567 ENDDO 568 ENDIF 569 ENDDO 570 ENDDO 606 DO k = nzb_u_inner(j,i) + 1, nzt 607 u(k,j,i) = u(k,j,i) + volume_flow_offset(1) 608 v(k,j,i) = v(k,j,i) + volume_flow_offset(2) 609 ENDDO 610 ENDDO 611 ENDDO 612 571 613 !$OMP END PARALLEL 572 614 … … 575 617 ! 576 618 ! Exchange of boundaries for the velocities 577 CALL exchange_horiz( u )578 CALL exchange_horiz( v )579 CALL exchange_horiz( w )619 CALL exchange_horiz( u, nbgp ) 620 CALL exchange_horiz( v, nbgp ) 621 CALL exchange_horiz( w, nbgp ) 580 622 581 623 ! … … 620 662 ENDDO 621 663 #endif 664 622 665 localsum = localsum + threadsum 623 666 !$OMP END PARALLEL … … 631 674 632 675 CALL cpu_log( log_point(8), 'pres', 'stop' ) 676 633 677 634 678
Note: See TracChangeset
for help on using the changeset viewer.