Changeset 978 for palm/trunk/SOURCE/boundary_conds.f90
- Timestamp:
- Aug 9, 2012 8:28:32 AM (12 years ago)
- Location:
- palm/trunk
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk
- Property svn:mergeinfo changed
/palm/branches/fricke (added) merged: 942-944,967-968,971-972,977
- Property svn:mergeinfo changed
-
palm/trunk/SOURCE
- Property svn:mergeinfo changed
/palm/branches/fricke/SOURCE (added) merged: 967-968,971-972,977
- Property svn:mergeinfo changed
-
palm/trunk/SOURCE/boundary_conds.f90
r876 r978 4 4 ! Current revisions: 5 5 ! ----------------- 6 ! 6 ! Neumann boudnary conditions are added at the inflow boundary for the SGS-TKE. 7 ! Outflow boundary conditions for the velocity components can be set to Neumann 8 ! conditions or to radiation conditions with a horizontal averaged phase 9 ! velocity. 7 10 ! 8 11 ! … … 190 193 q_p(nzt+1,:,:) = q_p(nzt,:,:) + bc_q_t_val * dzu(nzt+1) 191 194 192 193 195 ENDIF 194 196 … … 198 200 !-- Since in prognostic_equations (cache optimized version) these levels are 199 201 !-- handled as a prognostic level, boundary values have to be restored here. 202 !-- For the SGS-TKE, Neumann boundary conditions are used at the inflow. 200 203 IF ( inflow_s ) THEN 201 204 v_p(:,nys,:) = v_p(:,nys-1,:) 205 IF ( .NOT. constant_diffusion ) e_p(:,nys-1,:) = e_p(:,nys,:) 206 ELSEIF ( inflow_n ) THEN 207 IF ( .NOT. constant_diffusion ) e_p(:,nyn+1,:) = e_p(:,nyn,:) 202 208 ELSEIF ( inflow_l ) THEN 203 209 u_p(:,:,nxl) = u_p(:,:,nxl-1) 210 IF ( .NOT. constant_diffusion ) e_p(:,:,nxl-1) = e_p(:,:,nxl) 211 ELSEIF ( inflow_r ) THEN 212 IF ( .NOT. constant_diffusion ) e_p(:,:,nxr+1) = e_p(:,:,nxr) 204 213 ENDIF 205 214 … … 227 236 228 237 ! 229 !-- Radiation boundary condition for the velocities at the respective outflow 238 !-- Neumann or Radiation boundary condition for the velocities at the 239 !-- respective outflow 230 240 IF ( outflow_s ) THEN 231 241 232 c_max = dy / dt_3d 233 234 DO i = nxlg, nxrg 242 IF ( bc_ns_dirneu ) THEN 243 u(:,-1,:) = u(:,0,:) 244 v(:,0,:) = v(:,1,:) 245 w(:,-1,:) = w(:,0,:) 246 ELSEIF ( bc_ns_dirrad ) THEN 247 248 c_max = dy / dt_3d 249 250 c_u_m_l = 0.0 251 c_v_m_l = 0.0 252 c_w_m_l = 0.0 253 254 c_u_m = 0.0 255 c_v_m = 0.0 256 c_w_m = 0.0 257 258 ! 259 !-- Calculate the phase speeds for u,v, and w, first local and then 260 !-- average parallel along the outflow boundary. 235 261 DO k = nzb+1, nzt+1 236 237 ! 238 !-- Calculate the phase speeds for u,v, and w. In case of using a 239 !-- Runge-Kutta scheme, do this for the first substep only and then 240 !-- reuse this values for the further substeps. 241 IF ( intermediate_timestep_count == 1 ) THEN 262 DO i = nxl, nxr 242 263 243 264 denom = u_m_s(k,0,i) - u_m_s(k,1,i) 244 265 245 266 IF ( denom /= 0.0 ) THEN 246 c_u(k,i) = -c_max * ( u(k,0,i) - u_m_s(k,0,i) ) / denom 267 c_u(k,i) = -c_max * ( u(k,0,i) - u_m_s(k,0,i) ) & 268 / ( denom * tsc(2) ) 247 269 IF ( c_u(k,i) < 0.0 ) THEN 248 270 c_u(k,i) = 0.0 … … 257 279 258 280 IF ( denom /= 0.0 ) THEN 259 c_v(k,i) = -c_max * ( v(k,1,i) - v_m_s(k,1,i) ) / denom 281 c_v(k,i) = -c_max * ( v(k,1,i) - v_m_s(k,1,i) ) & 282 / ( denom * tsc(2) ) 260 283 IF ( c_v(k,i) < 0.0 ) THEN 261 284 c_v(k,i) = 0.0 … … 270 293 271 294 IF ( denom /= 0.0 ) THEN 272 c_w(k,i) = -c_max * ( w(k,0,i) - w_m_s(k,0,i) ) / denom 295 c_w(k,i) = -c_max * ( w(k,0,i) - w_m_s(k,0,i) ) & 296 / ( denom * tsc(2) ) 273 297 IF ( c_w(k,i) < 0.0 ) THEN 274 298 c_w(k,i) = 0.0 … … 280 304 ENDIF 281 305 282 ! 283 !-- Save old timelevels for the next timestep 284 u_m_s(k,:,i) = u(k,0:1,i) 285 v_m_s(k,:,i) = v(k,1:2,i) 286 w_m_s(k,:,i) = w(k,0:1,i) 287 288 ENDIF 289 290 ! 291 !-- Calculate the new velocities 292 u_p(k,-1,i) = u(k,-1,i) - dt_3d * tsc(2) * c_u(k,i) * & 306 c_u_m_l(k) = c_u_m_l(k) + c_u(k,i) 307 c_v_m_l(k) = c_v_m_l(k) + c_v(k,i) 308 c_w_m_l(k) = c_w_m_l(k) + c_w(k,i) 309 310 ENDDO 311 ENDDO 312 313 #if defined( __parallel ) 314 IF ( collective_wait ) CALL MPI_BARRIER( comm1dx, ierr ) 315 CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, & 316 MPI_SUM, comm1dx, ierr ) 317 IF ( collective_wait ) CALL MPI_BARRIER( comm1dx, ierr ) 318 CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, & 319 MPI_SUM, comm1dx, ierr ) 320 IF ( collective_wait ) CALL MPI_BARRIER( comm1dx, ierr ) 321 CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, & 322 MPI_SUM, comm1dx, ierr ) 323 #else 324 c_u_m = c_u_m_l 325 c_v_m = c_v_m_l 326 c_w_m = c_w_m_l 327 #endif 328 329 c_u_m = c_u_m / (nx+1) 330 c_v_m = c_v_m / (nx+1) 331 c_w_m = c_w_m / (nx+1) 332 333 ! 334 !-- Save old timelevels for the next timestep 335 IF ( intermediate_timestep_count == 1 ) THEN 336 u_m_s(:,:,:) = u(:,0:1,:) 337 v_m_s(:,:,:) = v(:,1:2,:) 338 w_m_s(:,:,:) = w(:,0:1,:) 339 ENDIF 340 341 ! 342 !-- Calculate the new velocities 343 DO k = nzb+1, nzt+1 344 DO i = nxlg, nxrg 345 u_p(k,-1,i) = u(k,-1,i) - dt_3d * tsc(2) * c_u_m(k) * & 293 346 ( u(k,-1,i) - u(k,0,i) ) * ddy 294 347 295 v_p(k,0,i) = v(k,0,i) - dt_3d * tsc(2) * c_v(k,i) *&348 v_p(k,0,i) = v(k,0,i) - dt_3d * tsc(2) * c_v_m(k) * & 296 349 ( v(k,0,i) - v(k,1,i) ) * ddy 297 350 298 w_p(k,-1,i) = w(k,-1,i) - dt_3d * tsc(2) * c_w(k,i) *&351 w_p(k,-1,i) = w(k,-1,i) - dt_3d * tsc(2) * c_w_m(k) * & 299 352 ( w(k,-1,i) - w(k,0,i) ) * ddy 300 301 ENDDO 302 ENDDO 303 304 ! 305 !-- Bottom boundary at the outflow 306 IF ( ibc_uv_b == 0 ) THEN 307 u_p(nzb,-1,:) = 0.0 308 v_p(nzb,0,:) = 0.0 309 ELSE 310 u_p(nzb,-1,:) = u_p(nzb+1,-1,:) 311 v_p(nzb,0,:) = v_p(nzb+1,0,:) 312 ENDIF 313 w_p(nzb,-1,:) = 0.0 314 315 ! 316 !-- Top boundary at the outflow 317 IF ( ibc_uv_t == 0 ) THEN 318 u_p(nzt+1,-1,:) = u_init(nzt+1) 319 v_p(nzt+1,0,:) = v_init(nzt+1) 320 ELSE 321 u_p(nzt+1,-1,:) = u(nzt,-1,:) 322 v_p(nzt+1,0,:) = v(nzt,0,:) 323 ENDIF 324 w_p(nzt:nzt+1,-1,:) = 0.0 353 ENDDO 354 ENDDO 355 356 ! 357 !-- Bottom boundary at the outflow 358 IF ( ibc_uv_b == 0 ) THEN 359 u_p(nzb,-1,:) = 0.0 360 v_p(nzb,0,:) = 0.0 361 ELSE 362 u_p(nzb,-1,:) = u_p(nzb+1,-1,:) 363 v_p(nzb,0,:) = v_p(nzb+1,0,:) 364 ENDIF 365 w_p(nzb,-1,:) = 0.0 366 367 ! 368 !-- Top boundary at the outflow 369 IF ( ibc_uv_t == 0 ) THEN 370 u_p(nzt+1,-1,:) = u_init(nzt+1) 371 v_p(nzt+1,0,:) = v_init(nzt+1) 372 ELSE 373 u_p(nzt+1,-1,:) = u(nzt,-1,:) 374 v_p(nzt+1,0,:) = v(nzt,0,:) 375 ENDIF 376 w_p(nzt:nzt+1,-1,:) = 0.0 377 378 ENDIF 325 379 326 380 ENDIF … … 328 382 IF ( outflow_n ) THEN 329 383 330 c_max = dy / dt_3d 331 332 DO i = nxlg, nxrg 384 IF ( bc_ns_neudir ) THEN 385 u(:,ny+1,:) = u(:,ny,:) 386 v(:,ny+1,:) = v(:,ny,:) 387 w(:,ny+1,:) = w(:,ny,:) 388 ELSEIF ( bc_ns_dirrad ) THEN 389 390 c_max = dy / dt_3d 391 392 c_u_m_l = 0.0 393 c_v_m_l = 0.0 394 c_w_m_l = 0.0 395 396 c_u_m = 0.0 397 c_v_m = 0.0 398 c_w_m = 0.0 399 400 ! 401 !-- Calculate the phase speeds for u,v, and w, first local and then 402 !-- average parallel along the outflow boundary. 333 403 DO k = nzb+1, nzt+1 334 335 ! 336 !-- Calculate the phase speeds for u,v, and w. In case of using a 337 !-- Runge-Kutta scheme, do this for the first substep only and then 338 !-- reuse this values for the further substeps. 339 IF ( intermediate_timestep_count == 1 ) THEN 404 DO i = nxl, nxr 340 405 341 406 denom = u_m_n(k,ny,i) - u_m_n(k,ny-1,i) 342 407 343 408 IF ( denom /= 0.0 ) THEN 344 c_u(k,i) = -c_max * ( u(k,ny,i) - u_m_n(k,ny,i) ) / denom 409 c_u(k,i) = -c_max * ( u(k,ny,i) - u_m_n(k,ny,i) ) & 410 / ( denom * tsc(2) ) 345 411 IF ( c_u(k,i) < 0.0 ) THEN 346 412 c_u(k,i) = 0.0 … … 355 421 356 422 IF ( denom /= 0.0 ) THEN 357 c_v(k,i) = -c_max * ( v(k,ny,i) - v_m_n(k,ny,i) ) / denom 423 c_v(k,i) = -c_max * ( v(k,ny,i) - v_m_n(k,ny,i) ) & 424 / ( denom * tsc(2) ) 358 425 IF ( c_v(k,i) < 0.0 ) THEN 359 426 c_v(k,i) = 0.0 … … 368 435 369 436 IF ( denom /= 0.0 ) THEN 370 c_w(k,i) = -c_max * ( w(k,ny,i) - w_m_n(k,ny,i) ) / denom 437 c_w(k,i) = -c_max * ( w(k,ny,i) - w_m_n(k,ny,i) ) & 438 / ( denom * tsc(2) ) 371 439 IF ( c_w(k,i) < 0.0 ) THEN 372 440 c_w(k,i) = 0.0 … … 378 446 ENDIF 379 447 380 ! 381 !-- Swap timelevels for the next timestep 382 u_m_n(k,:,i) = u(k,ny-1:ny,i) 383 v_m_n(k,:,i) = v(k,ny-1:ny,i) 384 w_m_n(k,:,i) = w(k,ny-1:ny,i) 385 386 ENDIF 387 388 ! 389 !-- Calculate the new velocities 390 u_p(k,ny+1,i) = u(k,ny+1,i) - dt_3d * tsc(2) * c_u(k,i) * & 391 ( u(k,ny+1,i) - u(k,ny,i) ) * ddy 392 393 v_p(k,ny+1,i) = v(k,ny+1,i) - dt_3d * tsc(2) * c_v(k,i) * & 394 ( v(k,ny+1,i) - v(k,ny,i) ) * ddy 395 396 w_p(k,ny+1,i) = w(k,ny+1,i) - dt_3d * tsc(2) * c_w(k,i) * & 397 ( w(k,ny+1,i) - w(k,ny,i) ) * ddy 398 399 ENDDO 400 ENDDO 401 402 ! 403 !-- Bottom boundary at the outflow 404 IF ( ibc_uv_b == 0 ) THEN 405 u_p(nzb,ny+1,:) = 0.0 406 v_p(nzb,ny+1,:) = 0.0 407 ELSE 408 u_p(nzb,ny+1,:) = u_p(nzb+1,ny+1,:) 409 v_p(nzb,ny+1,:) = v_p(nzb+1,ny+1,:) 410 ENDIF 411 w_p(nzb,ny+1,:) = 0.0 412 413 ! 414 !-- Top boundary at the outflow 415 IF ( ibc_uv_t == 0 ) THEN 416 u_p(nzt+1,ny+1,:) = u_init(nzt+1) 417 v_p(nzt+1,ny+1,:) = v_init(nzt+1) 418 ELSE 419 u_p(nzt+1,ny+1,:) = u_p(nzt,nyn+1,:) 420 v_p(nzt+1,ny+1,:) = v_p(nzt,nyn+1,:) 421 ENDIF 422 w_p(nzt:nzt+1,ny+1,:) = 0.0 448 c_u_m_l(k) = c_u_m_l(k) + c_u(k,i) 449 c_v_m_l(k) = c_v_m_l(k) + c_v(k,i) 450 c_w_m_l(k) = c_w_m_l(k) + c_w(k,i) 451 452 ENDDO 453 ENDDO 454 455 #if defined( __parallel ) 456 IF ( collective_wait ) CALL MPI_BARRIER( comm1dx, ierr ) 457 CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, & 458 MPI_SUM, comm1dx, ierr ) 459 IF ( collective_wait ) CALL MPI_BARRIER( comm1dx, ierr ) 460 CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, & 461 MPI_SUM, comm1dx, ierr ) 462 IF ( collective_wait ) CALL MPI_BARRIER( comm1dx, ierr ) 463 CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, & 464 MPI_SUM, comm1dx, ierr ) 465 #else 466 c_u_m = c_u_m_l 467 c_v_m = c_v_m_l 468 c_w_m = c_w_m_l 469 #endif 470 471 c_u_m = c_u_m / (nx+1) 472 c_v_m = c_v_m / (nx+1) 473 c_w_m = c_w_m / (nx+1) 474 475 ! 476 !-- Save old timelevels for the next timestep 477 IF ( intermediate_timestep_count == 1 ) THEN 478 u_m_n(:,:,:) = u(:,ny-1:ny,:) 479 v_m_n(:,:,:) = v(:,ny-1:ny,:) 480 w_m_n(:,:,:) = w(:,ny-1:ny,:) 481 ENDIF 482 483 ! 484 !-- Calculate the new velocities 485 DO k = nzb+1, nzt+1 486 DO i = nxlg, nxrg 487 u_p(k,ny+1,i) = u(k,ny+1,i) - dt_3d * tsc(2) * c_u_m(k) * & 488 ( u(k,ny+1,i) - u(k,ny,i) ) * ddy 489 490 v_p(k,ny+1,i) = v(k,ny+1,i) - dt_3d * tsc(2) * c_v_m(k) * & 491 ( v(k,ny+1,i) - v(k,ny,i) ) * ddy 492 493 w_p(k,ny+1,i) = w(k,ny+1,i) - dt_3d * tsc(2) * c_w_m(k) * & 494 ( w(k,ny+1,i) - w(k,ny,i) ) * ddy 495 ENDDO 496 ENDDO 497 498 ! 499 !-- Bottom boundary at the outflow 500 IF ( ibc_uv_b == 0 ) THEN 501 u_p(nzb,ny+1,:) = 0.0 502 v_p(nzb,ny+1,:) = 0.0 503 ELSE 504 u_p(nzb,ny+1,:) = u_p(nzb+1,ny+1,:) 505 v_p(nzb,ny+1,:) = v_p(nzb+1,ny+1,:) 506 ENDIF 507 w_p(nzb,ny+1,:) = 0.0 508 509 ! 510 !-- Top boundary at the outflow 511 IF ( ibc_uv_t == 0 ) THEN 512 u_p(nzt+1,ny+1,:) = u_init(nzt+1) 513 v_p(nzt+1,ny+1,:) = v_init(nzt+1) 514 ELSE 515 u_p(nzt+1,ny+1,:) = u_p(nzt,nyn+1,:) 516 v_p(nzt+1,ny+1,:) = v_p(nzt,nyn+1,:) 517 ENDIF 518 w_p(nzt:nzt+1,ny+1,:) = 0.0 519 520 ENDIF 423 521 424 522 ENDIF … … 426 524 IF ( outflow_l ) THEN 427 525 428 c_max = dx / dt_3d 429 430 DO j = nysg, nyng 526 IF ( bc_lr_neudir ) THEN 527 u(:,:,-1) = u(:,:,0) 528 v(:,:,0) = v(:,:,1) 529 w(:,:,-1) = w(:,:,0) 530 ELSEIF ( bc_ns_dirrad ) THEN 531 532 c_max = dx / dt_3d 533 534 c_u_m_l = 0.0 535 c_v_m_l = 0.0 536 c_w_m_l = 0.0 537 538 c_u_m = 0.0 539 c_v_m = 0.0 540 c_w_m = 0.0 541 542 ! 543 !-- Calculate the phase speeds for u,v, and w, first local and then 544 !-- average parallel along the outflow boundary. 431 545 DO k = nzb+1, nzt+1 432 433 ! 434 !-- Calculate the phase speeds for u,v, and w. In case of using a 435 !-- Runge-Kutta scheme, do this for the first substep only and then 436 !-- reuse this values for the further substeps. 437 IF ( intermediate_timestep_count == 1 ) THEN 546 DO j = nys, nyn 438 547 439 548 denom = u_m_l(k,j,1) - u_m_l(k,j,2) 440 549 441 550 IF ( denom /= 0.0 ) THEN 442 c_u(k,j) = -c_max * ( u(k,j,1) - u_m_l(k,j,1) ) / denom 551 c_u(k,j) = -c_max * ( u(k,j,1) - u_m_l(k,j,1) ) & 552 / ( denom * tsc(2) ) 443 553 IF ( c_u(k,j) < 0.0 ) THEN 444 554 c_u(k,j) = 0.0 … … 453 563 454 564 IF ( denom /= 0.0 ) THEN 455 c_v(k,j) = -c_max * ( v(k,j,0) - v_m_l(k,j,0) ) / denom 565 c_v(k,j) = -c_max * ( v(k,j,0) - v_m_l(k,j,0) ) & 566 / ( denom * tsc(2) ) 456 567 IF ( c_v(k,j) < 0.0 ) THEN 457 568 c_v(k,j) = 0.0 … … 466 577 467 578 IF ( denom /= 0.0 ) THEN 468 c_w(k,j) = -c_max * ( w(k,j,0) - w_m_l(k,j,0) ) / denom 579 c_w(k,j) = -c_max * ( w(k,j,0) - w_m_l(k,j,0) ) & 580 / ( denom * tsc(2) ) 469 581 IF ( c_w(k,j) < 0.0 ) THEN 470 582 c_w(k,j) = 0.0 … … 476 588 ENDIF 477 589 478 ! 479 !-- Swap timelevels for the next timestep 480 u_m_l(k,j,:) = u(k,j,1:2) 481 v_m_l(k,j,:) = v(k,j,0:1) 482 w_m_l(k,j,:) = w(k,j,0:1) 483 484 ENDIF 485 486 ! 487 !-- Calculate the new velocities 488 u_p(k,j,0) = u(k,j,0) - dt_3d * tsc(2) * c_u(k,j) * & 590 c_u_m_l(k) = c_u_m_l(k) + c_u(k,j) 591 c_v_m_l(k) = c_v_m_l(k) + c_v(k,j) 592 c_w_m_l(k) = c_w_m_l(k) + c_w(k,j) 593 594 ENDDO 595 ENDDO 596 597 #if defined( __parallel ) 598 IF ( collective_wait ) CALL MPI_BARRIER( comm1dy, ierr ) 599 CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, & 600 MPI_SUM, comm1dy, ierr ) 601 IF ( collective_wait ) CALL MPI_BARRIER( comm1dy, ierr ) 602 CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, & 603 MPI_SUM, comm1dy, ierr ) 604 IF ( collective_wait ) CALL MPI_BARRIER( comm1dy, ierr ) 605 CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, & 606 MPI_SUM, comm1dy, ierr ) 607 #else 608 c_u_m = c_u_m_l 609 c_v_m = c_v_m_l 610 c_w_m = c_w_m_l 611 #endif 612 613 c_u_m = c_u_m / (ny+1) 614 c_v_m = c_v_m / (ny+1) 615 c_w_m = c_w_m / (ny+1) 616 617 ! 618 !-- Save old timelevels for the next timestep 619 IF ( intermediate_timestep_count == 1 ) THEN 620 u_m_l(:,:,:) = u(:,:,1:2) 621 v_m_l(:,:,:) = v(:,:,0:1) 622 w_m_l(:,:,:) = w(:,:,0:1) 623 ENDIF 624 625 ! 626 !-- Calculate the new velocities 627 DO k = nzb+1, nzt+1 628 DO i = nxlg, nxrg 629 u_p(k,j,0) = u(k,j,0) - dt_3d * tsc(2) * c_u_m(k) * & 489 630 ( u(k,j,0) - u(k,j,1) ) * ddx 490 631 491 v_p(k,j,-1) = v(k,j,-1) - dt_3d * tsc(2) * c_v(k,j) *&632 v_p(k,j,-1) = v(k,j,-1) - dt_3d * tsc(2) * c_v_m(k) * & 492 633 ( v(k,j,-1) - v(k,j,0) ) * ddx 493 634 494 w_p(k,j,-1) = w(k,j,-1) - dt_3d * tsc(2) * c_w(k,j) *&635 w_p(k,j,-1) = w(k,j,-1) - dt_3d * tsc(2) * c_w_m(k) * & 495 636 ( w(k,j,-1) - w(k,j,0) ) * ddx 496 497 ENDDO 498 ENDDO 499 500 ! 501 !-- Bottom boundary at the outflow 502 IF ( ibc_uv_b == 0 ) THEN 503 u_p(nzb,:,0) = 0.0 504 v_p(nzb,:,-1) = 0.0 505 ELSE 506 u_p(nzb,:,0) = u_p(nzb+1,:,0) 507 v_p(nzb,:,-1) = v_p(nzb+1,:,-1) 508 ENDIF 509 w_p(nzb,:,-1) = 0.0 510 511 ! 512 !-- Top boundary at the outflow 513 IF ( ibc_uv_t == 0 ) THEN 514 u_p(nzt+1,:,-1) = u_init(nzt+1) 515 v_p(nzt+1,:,-1) = v_init(nzt+1) 516 ELSE 517 u_p(nzt+1,:,-1) = u_p(nzt,:,-1) 518 v_p(nzt+1,:,-1) = v_p(nzt,:,-1) 519 ENDIF 520 w_p(nzt:nzt+1,:,-1) = 0.0 637 ENDDO 638 ENDDO 639 640 ! 641 !-- Bottom boundary at the outflow 642 IF ( ibc_uv_b == 0 ) THEN 643 u_p(nzb,:,0) = 0.0 644 v_p(nzb,:,-1) = 0.0 645 ELSE 646 u_p(nzb,:,0) = u_p(nzb+1,:,0) 647 v_p(nzb,:,-1) = v_p(nzb+1,:,-1) 648 ENDIF 649 w_p(nzb,:,-1) = 0.0 650 651 ! 652 !-- Top boundary at the outflow 653 IF ( ibc_uv_t == 0 ) THEN 654 u_p(nzt+1,:,-1) = u_init(nzt+1) 655 v_p(nzt+1,:,-1) = v_init(nzt+1) 656 ELSE 657 u_p(nzt+1,:,-1) = u_p(nzt,:,-1) 658 v_p(nzt+1,:,-1) = v_p(nzt,:,-1) 659 ENDIF 660 w_p(nzt:nzt+1,:,-1) = 0.0 661 662 ENDIF 521 663 522 664 ENDIF … … 524 666 IF ( outflow_r ) THEN 525 667 526 c_max = dx / dt_3d 527 528 DO j = nysg, nyng 668 IF ( bc_lr_dirneu ) THEN 669 u(:,:,nx+1) = u(:,:,nx) 670 v(:,:,nx+1) = v(:,:,nx) 671 w(:,:,nx+1) = w(:,:,nx) 672 ELSEIF ( bc_ns_dirrad ) THEN 673 674 c_max = dx / dt_3d 675 676 c_u_m_l = 0.0 677 c_v_m_l = 0.0 678 c_w_m_l = 0.0 679 680 c_u_m = 0.0 681 c_v_m = 0.0 682 c_w_m = 0.0 683 684 ! 685 !-- Calculate the phase speeds for u,v, and w, first local and then 686 !-- average parallel along the outflow boundary. 529 687 DO k = nzb+1, nzt+1 530 531 ! 532 !-- Calculate the phase speeds for u,v, and w. In case of using a 533 !-- Runge-Kutta scheme, do this for the first substep only and then 534 !-- reuse this values for the further substeps. 535 IF ( intermediate_timestep_count == 1 ) THEN 688 DO j = nys, nyn 536 689 537 690 denom = u_m_r(k,j,nx) - u_m_r(k,j,nx-1) 538 691 539 692 IF ( denom /= 0.0 ) THEN 540 c_u(k,j) = -c_max * ( u(k,j,nx) - u_m_r(k,j,nx) ) / denom 693 c_u(k,j) = -c_max * ( u(k,j,nx) - u_m_r(k,j,nx) ) & 694 / ( denom * tsc(2) ) 541 695 IF ( c_u(k,j) < 0.0 ) THEN 542 696 c_u(k,j) = 0.0 … … 551 705 552 706 IF ( denom /= 0.0 ) THEN 553 c_v(k,j) = -c_max * ( v(k,j,nx) - v_m_r(k,j,nx) ) / denom 707 c_v(k,j) = -c_max * ( v(k,j,nx) - v_m_r(k,j,nx) ) & 708 / ( denom * tsc(2) ) 554 709 IF ( c_v(k,j) < 0.0 ) THEN 555 710 c_v(k,j) = 0.0 … … 564 719 565 720 IF ( denom /= 0.0 ) THEN 566 c_w(k,j) = -c_max * ( w(k,j,nx) - w_m_r(k,j,nx) ) / denom 721 c_w(k,j) = -c_max * ( w(k,j,nx) - w_m_r(k,j,nx) ) & 722 / ( denom * tsc(2) ) 567 723 IF ( c_w(k,j) < 0.0 ) THEN 568 724 c_w(k,j) = 0.0 … … 574 730 ENDIF 575 731 576 ! 577 !-- Swap timelevels for the next timestep 578 u_m_r(k,j,:) = u(k,j,nx-1:nx) 579 v_m_r(k,j,:) = v(k,j,nx-1:nx) 580 w_m_r(k,j,:) = w(k,j,nx-1:nx) 581 582 ENDIF 583 584 ! 585 !-- Calculate the new velocities 586 u_p(k,j,nx+1) = u(k,j,nx+1) - dt_3d * tsc(2) * c_u(k,j) * & 587 ( u(k,j,nx+1) - u(k,j,nx) ) * ddx 588 589 v_p(k,j,nx+1) = v(k,j,nx+1) - dt_3d * tsc(2) * c_v(k,j) * & 590 ( v(k,j,nx+1) - v(k,j,nx) ) * ddx 591 592 w_p(k,j,nx+1) = w(k,j,nx+1) - dt_3d * tsc(2) * c_w(k,j) * & 593 ( w(k,j,nx+1) - w(k,j,nx) ) * ddx 594 595 ENDDO 596 ENDDO 597 598 ! 599 !-- Bottom boundary at the outflow 600 IF ( ibc_uv_b == 0 ) THEN 601 u_p(nzb,:,nx+1) = 0.0 602 v_p(nzb,:,nx+1) = 0.0 603 ELSE 604 u_p(nzb,:,nx+1) = u_p(nzb+1,:,nx+1) 605 v_p(nzb,:,nx+1) = v_p(nzb+1,:,nx+1) 606 ENDIF 607 w_p(nzb,:,nx+1) = 0.0 608 609 ! 610 !-- Top boundary at the outflow 611 IF ( ibc_uv_t == 0 ) THEN 612 u_p(nzt+1,:,nx+1) = u_init(nzt+1) 613 v_p(nzt+1,:,nx+1) = v_init(nzt+1) 614 ELSE 615 u_p(nzt+1,:,nx+1) = u_p(nzt,:,nx+1) 616 v_p(nzt+1,:,nx+1) = v_p(nzt,:,nx+1) 617 ENDIF 618 w(nzt:nzt+1,:,nx+1) = 0.0 732 c_u_m_l(k) = c_u_m_l(k) + c_u(k,j) 733 c_v_m_l(k) = c_v_m_l(k) + c_v(k,j) 734 c_w_m_l(k) = c_w_m_l(k) + c_w(k,j) 735 736 ENDDO 737 ENDDO 738 739 #if defined( __parallel ) 740 IF ( collective_wait ) CALL MPI_BARRIER( comm1dy, ierr ) 741 CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, & 742 MPI_SUM, comm1dy, ierr ) 743 IF ( collective_wait ) CALL MPI_BARRIER( comm1dy, ierr ) 744 CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, & 745 MPI_SUM, comm1dy, ierr ) 746 IF ( collective_wait ) CALL MPI_BARRIER( comm1dy, ierr ) 747 CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, & 748 MPI_SUM, comm1dy, ierr ) 749 #else 750 c_u_m = c_u_m_l 751 c_v_m = c_v_m_l 752 c_w_m = c_w_m_l 753 #endif 754 755 c_u_m = c_u_m / (ny+1) 756 c_v_m = c_v_m / (ny+1) 757 c_w_m = c_w_m / (ny+1) 758 759 ! 760 !-- Save old timelevels for the next timestep 761 IF ( intermediate_timestep_count == 1 ) THEN 762 u_m_r(:,:,:) = u(:,:,nx-1:nx) 763 v_m_r(:,:,:) = v(:,:,nx-1:nx) 764 w_m_r(:,:,:) = w(:,:,nx-1:nx) 765 ENDIF 766 767 ! 768 !-- Calculate the new velocities 769 DO k = nzb+1, nzt+1 770 DO i = nxlg, nxrg 771 u_p(k,j,nx+1) = u(k,j,nx+1) - dt_3d * tsc(2) * c_u_m(k) * & 772 ( u(k,j,nx+1) - u(k,j,nx) ) * ddx 773 774 v_p(k,j,nx+1) = v(k,j,nx+1) - dt_3d * tsc(2) * c_v_m(k) * & 775 ( v(k,j,nx+1) - v(k,j,nx) ) * ddx 776 777 w_p(k,j,nx+1) = w(k,j,nx+1) - dt_3d * tsc(2) * c_w_m(k) * & 778 ( w(k,j,nx+1) - w(k,j,nx) ) * ddx 779 ENDDO 780 ENDDO 781 782 ! 783 !-- Bottom boundary at the outflow 784 IF ( ibc_uv_b == 0 ) THEN 785 u_p(nzb,:,nx+1) = 0.0 786 v_p(nzb,:,nx+1) = 0.0 787 ELSE 788 u_p(nzb,:,nx+1) = u_p(nzb+1,:,nx+1) 789 v_p(nzb,:,nx+1) = v_p(nzb+1,:,nx+1) 790 ENDIF 791 w_p(nzb,:,nx+1) = 0.0 792 793 ! 794 !-- Top boundary at the outflow 795 IF ( ibc_uv_t == 0 ) THEN 796 u_p(nzt+1,:,nx+1) = u_init(nzt+1) 797 v_p(nzt+1,:,nx+1) = v_init(nzt+1) 798 ELSE 799 u_p(nzt+1,:,nx+1) = u_p(nzt,:,nx+1) 800 v_p(nzt+1,:,nx+1) = v_p(nzt,:,nx+1) 801 ENDIF 802 w(nzt:nzt+1,:,nx+1) = 0.0 803 804 ENDIF 619 805 620 806 ENDIF 621 807 622 623 808 END SUBROUTINE boundary_conds
Note: See TracChangeset
for help on using the changeset viewer.
