Changeset 4268 for palm/trunk/SOURCE/boundary_conds.f90
- Timestamp:
- Oct 17, 2019 11:29:38 AM (5 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/boundary_conds.f90 (modified) (23 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/boundary_conds.f90
r4182 r4268 25 25 ! ----------------- 26 26 ! $Id$ 27 ! Removing bulk cloud variables to respective module 28 ! 29 ! 4182 2019-08-22 15:20:23Z scharf 27 30 ! Corrected "Former revisions" section 28 31 ! … … 64 67 USE arrays_3d, & 65 68 ONLY: c_u, c_u_m, c_u_m_l, c_v, c_v_m, c_v_m_l, c_w, c_w_m, c_w_m_l, & 66 dzu, nc_p, nr_p, pt, pt_init, pt_p, q,&67 q_p, qc_p, qr_p, s, s_p, sa, sa_p, u, u_init, u_m_l, u_m_n,&69 dzu, pt, pt_init, pt_p, q, & 70 q_p, s, s_p, sa, sa_p, u, u_init, u_m_l, u_m_n, & 68 71 u_m_r, u_m_s, u_p, v, v_init, v_m_l, v_m_n, v_m_r, v_m_s, v_p, & 69 72 w, w_p, w_m_l, w_m_n, w_m_r, w_m_s 70 71 USE bulk_cloud_model_mod, &72 ONLY: bulk_cloud_model, microphysics_morrison, microphysics_seifert73 73 74 74 USE chemistry_model_mod, & … … 102 102 103 103 USE salsa_mod, & 104 ONLY: salsa_boundary_conds 104 ONLY: salsa_boundary_conds 105 105 106 106 USE surface_mod, & … … 136 136 !$ACC PRESENT(bc_h, w_p) 137 137 DO m = 1, bc_h(l)%ns 138 i = bc_h(l)%i(m) 138 i = bc_h(l)%i(m) 139 139 j = bc_h(l)%j(m) 140 140 k = bc_h(l)%k(m) … … 174 174 !$OMP PARALLEL DO PRIVATE( i, j, k ) 175 175 DO m = 1, bc_h(l)%ns 176 i = bc_h(l)%i(m) 176 i = bc_h(l)%i(m) 177 177 j = bc_h(l)%j(m) 178 178 k = bc_h(l)%k(m) … … 223 223 !$OMP PARALLEL DO PRIVATE( i, j, k ) 224 224 DO m = 1, bc_h(l)%ns 225 i = bc_h(l)%i(m) 225 i = bc_h(l)%i(m) 226 226 j = bc_h(l)%j(m) 227 227 k = bc_h(l)%k(m) … … 253 253 !$OMP PARALLEL DO PRIVATE( i, j, k ) 254 254 DO m = 1, bc_h(l)%ns 255 i = bc_h(l)%i(m) 255 i = bc_h(l)%i(m) 256 256 j = bc_h(l)%j(m) 257 257 k = bc_h(l)%k(m) … … 259 259 ENDDO 260 260 ENDDO 261 261 262 262 ELSE 263 263 264 264 DO l = 0, 1 265 265 !$OMP PARALLEL DO PRIVATE( i, j, k ) 266 266 DO m = 1, bc_h(l)%ns 267 i = bc_h(l)%i(m) 267 i = bc_h(l)%i(m) 268 268 j = bc_h(l)%j(m) 269 269 k = bc_h(l)%k(m) … … 279 279 q_p(nzt+1,:,:) = q_p(nzt,:,:) + bc_q_t_val * dzu(nzt+1) 280 280 ENDIF 281 282 IF ( bulk_cloud_model .AND. microphysics_morrison ) THEN283 !284 !-- Surface conditions cloud water (Dirichlet)285 !-- Run loop over all non-natural and natural walls. Note, in wall-datatype286 !-- the k coordinate belongs to the atmospheric grid point, therefore, set287 !-- qr_p and nr_p at upward (k-1) and downward-facing (k+1) walls288 DO l = 0, 1289 !$OMP PARALLEL DO PRIVATE( i, j, k )290 DO m = 1, bc_h(l)%ns291 i = bc_h(l)%i(m)292 j = bc_h(l)%j(m)293 k = bc_h(l)%k(m)294 qc_p(k+bc_h(l)%koff,j,i) = 0.0_wp295 nc_p(k+bc_h(l)%koff,j,i) = 0.0_wp296 ENDDO297 ENDDO298 !299 !-- Top boundary condition for cloud water (Dirichlet)300 qc_p(nzt+1,:,:) = 0.0_wp301 nc_p(nzt+1,:,:) = 0.0_wp302 303 ENDIF304 305 IF ( bulk_cloud_model .AND. microphysics_seifert ) THEN306 !307 !-- Surface conditions rain water (Dirichlet)308 !-- Run loop over all non-natural and natural walls. Note, in wall-datatype309 !-- the k coordinate belongs to the atmospheric grid point, therefore, set310 !-- qr_p and nr_p at upward (k-1) and downward-facing (k+1) walls311 DO l = 0, 1312 !$OMP PARALLEL DO PRIVATE( i, j, k )313 DO m = 1, bc_h(l)%ns314 i = bc_h(l)%i(m)315 j = bc_h(l)%j(m)316 k = bc_h(l)%k(m)317 qr_p(k+bc_h(l)%koff,j,i) = 0.0_wp318 nr_p(k+bc_h(l)%koff,j,i) = 0.0_wp319 ENDDO320 ENDDO321 !322 !-- Top boundary condition for rain water (Dirichlet)323 qr_p(nzt+1,:,:) = 0.0_wp324 nr_p(nzt+1,:,:) = 0.0_wp325 326 ENDIF327 281 ENDIF 328 282 ! … … 336 290 !-- s_p at k-1 337 291 IF ( ibc_s_b == 0 ) THEN 338 292 339 293 DO l = 0, 1 340 294 !$OMP PARALLEL DO PRIVATE( i, j, k ) 341 295 DO m = 1, bc_h(l)%ns 342 i = bc_h(l)%i(m) 296 i = bc_h(l)%i(m) 343 297 j = bc_h(l)%j(m) 344 298 k = bc_h(l)%k(m) … … 346 300 ENDDO 347 301 ENDDO 348 302 349 303 ELSE 350 304 351 305 DO l = 0, 1 352 306 !$OMP PARALLEL DO PRIVATE( i, j, k ) 353 307 DO m = 1, bc_h(l)%ns 354 i = bc_h(l)%i(m) 308 i = bc_h(l)%i(m) 355 309 j = bc_h(l)%j(m) 356 310 k = bc_h(l)%k(m) … … 411 365 IF ( humidity ) THEN 412 366 q_p(:,nys-1,:) = q_p(:,nys,:) 413 IF ( bulk_cloud_model .AND. microphysics_morrison ) THEN414 qc_p(:,nys-1,:) = qc_p(:,nys,:)415 nc_p(:,nys-1,:) = nc_p(:,nys,:)416 ENDIF417 IF ( bulk_cloud_model .AND. microphysics_seifert ) THEN418 qr_p(:,nys-1,:) = qr_p(:,nys,:)419 nr_p(:,nys-1,:) = nr_p(:,nys,:)420 ENDIF421 367 ENDIF 422 368 IF ( passive_scalar ) s_p(:,nys-1,:) = s_p(:,nys,:) … … 425 371 IF ( humidity ) THEN 426 372 q_p(:,nyn+1,:) = q_p(:,nyn,:) 427 IF ( bulk_cloud_model .AND. microphysics_morrison ) THEN428 qc_p(:,nyn+1,:) = qc_p(:,nyn,:)429 nc_p(:,nyn+1,:) = nc_p(:,nyn,:)430 ENDIF431 IF ( bulk_cloud_model .AND. microphysics_seifert ) THEN432 qr_p(:,nyn+1,:) = qr_p(:,nyn,:)433 nr_p(:,nyn+1,:) = nr_p(:,nyn,:)434 ENDIF435 373 ENDIF 436 374 IF ( passive_scalar ) s_p(:,nyn+1,:) = s_p(:,nyn,:) … … 439 377 IF ( humidity ) THEN 440 378 q_p(:,:,nxl-1) = q_p(:,:,nxl) 441 IF ( bulk_cloud_model .AND. microphysics_morrison ) THEN442 qc_p(:,:,nxl-1) = qc_p(:,:,nxl)443 nc_p(:,:,nxl-1) = nc_p(:,:,nxl)444 ENDIF445 IF ( bulk_cloud_model .AND. microphysics_seifert ) THEN446 qr_p(:,:,nxl-1) = qr_p(:,:,nxl)447 nr_p(:,:,nxl-1) = nr_p(:,:,nxl)448 ENDIF449 379 ENDIF 450 380 IF ( passive_scalar ) s_p(:,:,nxl-1) = s_p(:,:,nxl) … … 453 383 IF ( humidity ) THEN 454 384 q_p(:,:,nxr+1) = q_p(:,:,nxr) 455 IF ( bulk_cloud_model .AND. microphysics_morrison ) THEN456 qc_p(:,:,nxr+1) = qc_p(:,:,nxr)457 nc_p(:,:,nxr+1) = nc_p(:,:,nxr)458 ENDIF459 IF ( bulk_cloud_model .AND. microphysics_seifert ) THEN460 qr_p(:,:,nxr+1) = qr_p(:,:,nxr)461 nr_p(:,:,nxr+1) = nr_p(:,:,nxr)462 ENDIF463 385 ENDIF 464 386 IF ( passive_scalar ) s_p(:,:,nxr+1) = s_p(:,:,nxr) … … 479 401 u_p(:,-1,:) = u(:,0,:) 480 402 v_p(:,0,:) = v(:,1,:) 481 w_p(:,-1,:) = w(:,0,:) 403 w_p(:,-1,:) = w(:,0,:) 482 404 ELSEIF ( .NOT. use_cmax ) THEN 483 405 … … 528 450 IF ( denom /= 0.0_wp ) THEN 529 451 c_w(k,i) = -c_max * ( w(k,0,i) - w_m_s(k,0,i) ) / ( denom * tsc(2) ) 452 IF ( c_w(k,i) < 0.0_wp ) THEN 453 c_w(k,i) = 0.0_wp 454 ELSEIF ( c_w(k,i) > c_max ) THEN 455 c_w(k,i) = c_max 456 ENDIF 457 ELSE 458 c_w(k,i) = c_max 459 ENDIF 460 461 c_u_m_l(k) = c_u_m_l(k) + c_u(k,i) 462 c_v_m_l(k) = c_v_m_l(k) + c_v(k,i) 463 c_w_m_l(k) = c_w_m_l(k) + c_w(k,i) 464 465 ENDDO 466 ENDDO 467 468 #if defined( __parallel ) 469 IF ( collective_wait ) CALL MPI_BARRIER( comm1dx, ierr ) 470 CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, & 471 MPI_SUM, comm1dx, ierr ) 472 IF ( collective_wait ) CALL MPI_BARRIER( comm1dx, ierr ) 473 CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, & 474 MPI_SUM, comm1dx, ierr ) 475 IF ( collective_wait ) CALL MPI_BARRIER( comm1dx, ierr ) 476 CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, & 477 MPI_SUM, comm1dx, ierr ) 478 #else 479 c_u_m = c_u_m_l 480 c_v_m = c_v_m_l 481 c_w_m = c_w_m_l 482 #endif 483 484 c_u_m = c_u_m / (nx+1) 485 c_v_m = c_v_m / (nx+1) 486 c_w_m = c_w_m / (nx+1) 487 488 ! 489 !-- Save old timelevels for the next timestep 490 IF ( intermediate_timestep_count == 1 ) THEN 491 u_m_s(:,:,:) = u(:,0:1,:) 492 v_m_s(:,:,:) = v(:,1:2,:) 493 w_m_s(:,:,:) = w(:,0:1,:) 494 ENDIF 495 496 ! 497 !-- Calculate the new velocities 498 DO k = nzb+1, nzt+1 499 DO i = nxlg, nxrg 500 u_p(k,-1,i) = u(k,-1,i) - dt_3d * tsc(2) * c_u_m(k) * & 501 ( u(k,-1,i) - u(k,0,i) ) * ddy 502 503 v_p(k,0,i) = v(k,0,i) - dt_3d * tsc(2) * c_v_m(k) * & 504 ( v(k,0,i) - v(k,1,i) ) * ddy 505 506 w_p(k,-1,i) = w(k,-1,i) - dt_3d * tsc(2) * c_w_m(k) * & 507 ( w(k,-1,i) - w(k,0,i) ) * ddy 508 ENDDO 509 ENDDO 510 511 ! 512 !-- Bottom boundary at the outflow 513 IF ( ibc_uv_b == 0 ) THEN 514 u_p(nzb,-1,:) = 0.0_wp 515 v_p(nzb,0,:) = 0.0_wp 516 ELSE 517 u_p(nzb,-1,:) = u_p(nzb+1,-1,:) 518 v_p(nzb,0,:) = v_p(nzb+1,0,:) 519 ENDIF 520 w_p(nzb,-1,:) = 0.0_wp 521 522 ! 523 !-- Top boundary at the outflow 524 IF ( ibc_uv_t == 0 ) THEN 525 u_p(nzt+1,-1,:) = u_init(nzt+1) 526 v_p(nzt+1,0,:) = v_init(nzt+1) 527 ELSE 528 u_p(nzt+1,-1,:) = u_p(nzt,-1,:) 529 v_p(nzt+1,0,:) = v_p(nzt,0,:) 530 ENDIF 531 w_p(nzt:nzt+1,-1,:) = 0.0_wp 532 533 ENDIF 534 535 ENDIF 536 537 IF ( bc_radiation_n ) THEN 538 539 IF ( use_cmax ) THEN 540 u_p(:,ny+1,:) = u(:,ny,:) 541 v_p(:,ny+1,:) = v(:,ny,:) 542 w_p(:,ny+1,:) = w(:,ny,:) 543 ELSEIF ( .NOT. use_cmax ) THEN 544 545 c_max = dy / dt_3d 546 547 c_u_m_l = 0.0_wp 548 c_v_m_l = 0.0_wp 549 c_w_m_l = 0.0_wp 550 551 c_u_m = 0.0_wp 552 c_v_m = 0.0_wp 553 c_w_m = 0.0_wp 554 555 ! 556 !-- Calculate the phase speeds for u, v, and w, first local and then 557 !-- average along the outflow boundary. 558 DO k = nzb+1, nzt+1 559 DO i = nxl, nxr 560 561 denom = u_m_n(k,ny,i) - u_m_n(k,ny-1,i) 562 563 IF ( denom /= 0.0_wp ) THEN 564 c_u(k,i) = -c_max * ( u(k,ny,i) - u_m_n(k,ny,i) ) / ( denom * tsc(2) ) 565 IF ( c_u(k,i) < 0.0_wp ) THEN 566 c_u(k,i) = 0.0_wp 567 ELSEIF ( c_u(k,i) > c_max ) THEN 568 c_u(k,i) = c_max 569 ENDIF 570 ELSE 571 c_u(k,i) = c_max 572 ENDIF 573 574 denom = v_m_n(k,ny,i) - v_m_n(k,ny-1,i) 575 576 IF ( denom /= 0.0_wp ) THEN 577 c_v(k,i) = -c_max * ( v(k,ny,i) - v_m_n(k,ny,i) ) / ( denom * tsc(2) ) 578 IF ( c_v(k,i) < 0.0_wp ) THEN 579 c_v(k,i) = 0.0_wp 580 ELSEIF ( c_v(k,i) > c_max ) THEN 581 c_v(k,i) = c_max 582 ENDIF 583 ELSE 584 c_v(k,i) = c_max 585 ENDIF 586 587 denom = w_m_n(k,ny,i) - w_m_n(k,ny-1,i) 588 589 IF ( denom /= 0.0_wp ) THEN 590 c_w(k,i) = -c_max * ( w(k,ny,i) - w_m_n(k,ny,i) ) / ( denom * tsc(2) ) 530 591 IF ( c_w(k,i) < 0.0_wp ) THEN 531 592 c_w(k,i) = 0.0_wp … … 567 628 !-- Save old timelevels for the next timestep 568 629 IF ( intermediate_timestep_count == 1 ) THEN 569 u_m_s(:,:,:) = u(:,0:1,:)570 v_m_s(:,:,:) = v(:,1:2,:)571 w_m_s(:,:,:) = w(:,0:1,:)572 ENDIF573 574 !575 !-- Calculate the new velocities576 DO k = nzb+1, nzt+1577 DO i = nxlg, nxrg578 u_p(k,-1,i) = u(k,-1,i) - dt_3d * tsc(2) * c_u_m(k) * &579 ( u(k,-1,i) - u(k,0,i) ) * ddy580 581 v_p(k,0,i) = v(k,0,i) - dt_3d * tsc(2) * c_v_m(k) * &582 ( v(k,0,i) - v(k,1,i) ) * ddy583 584 w_p(k,-1,i) = w(k,-1,i) - dt_3d * tsc(2) * c_w_m(k) * &585 ( w(k,-1,i) - w(k,0,i) ) * ddy586 ENDDO587 ENDDO588 589 !590 !-- Bottom boundary at the outflow591 IF ( ibc_uv_b == 0 ) THEN592 u_p(nzb,-1,:) = 0.0_wp593 v_p(nzb,0,:) = 0.0_wp594 ELSE595 u_p(nzb,-1,:) = u_p(nzb+1,-1,:)596 v_p(nzb,0,:) = v_p(nzb+1,0,:)597 ENDIF598 w_p(nzb,-1,:) = 0.0_wp599 600 !601 !-- Top boundary at the outflow602 IF ( ibc_uv_t == 0 ) THEN603 u_p(nzt+1,-1,:) = u_init(nzt+1)604 v_p(nzt+1,0,:) = v_init(nzt+1)605 ELSE606 u_p(nzt+1,-1,:) = u_p(nzt,-1,:)607 v_p(nzt+1,0,:) = v_p(nzt,0,:)608 ENDIF609 w_p(nzt:nzt+1,-1,:) = 0.0_wp610 611 ENDIF612 613 ENDIF614 615 IF ( bc_radiation_n ) THEN616 617 IF ( use_cmax ) THEN618 u_p(:,ny+1,:) = u(:,ny,:)619 v_p(:,ny+1,:) = v(:,ny,:)620 w_p(:,ny+1,:) = w(:,ny,:)621 ELSEIF ( .NOT. use_cmax ) THEN622 623 c_max = dy / dt_3d624 625 c_u_m_l = 0.0_wp626 c_v_m_l = 0.0_wp627 c_w_m_l = 0.0_wp628 629 c_u_m = 0.0_wp630 c_v_m = 0.0_wp631 c_w_m = 0.0_wp632 633 !634 !-- Calculate the phase speeds for u, v, and w, first local and then635 !-- average along the outflow boundary.636 DO k = nzb+1, nzt+1637 DO i = nxl, nxr638 639 denom = u_m_n(k,ny,i) - u_m_n(k,ny-1,i)640 641 IF ( denom /= 0.0_wp ) THEN642 c_u(k,i) = -c_max * ( u(k,ny,i) - u_m_n(k,ny,i) ) / ( denom * tsc(2) )643 IF ( c_u(k,i) < 0.0_wp ) THEN644 c_u(k,i) = 0.0_wp645 ELSEIF ( c_u(k,i) > c_max ) THEN646 c_u(k,i) = c_max647 ENDIF648 ELSE649 c_u(k,i) = c_max650 ENDIF651 652 denom = v_m_n(k,ny,i) - v_m_n(k,ny-1,i)653 654 IF ( denom /= 0.0_wp ) THEN655 c_v(k,i) = -c_max * ( v(k,ny,i) - v_m_n(k,ny,i) ) / ( denom * tsc(2) )656 IF ( c_v(k,i) < 0.0_wp ) THEN657 c_v(k,i) = 0.0_wp658 ELSEIF ( c_v(k,i) > c_max ) THEN659 c_v(k,i) = c_max660 ENDIF661 ELSE662 c_v(k,i) = c_max663 ENDIF664 665 denom = w_m_n(k,ny,i) - w_m_n(k,ny-1,i)666 667 IF ( denom /= 0.0_wp ) THEN668 c_w(k,i) = -c_max * ( w(k,ny,i) - w_m_n(k,ny,i) ) / ( denom * tsc(2) )669 IF ( c_w(k,i) < 0.0_wp ) THEN670 c_w(k,i) = 0.0_wp671 ELSEIF ( c_w(k,i) > c_max ) THEN672 c_w(k,i) = c_max673 ENDIF674 ELSE675 c_w(k,i) = c_max676 ENDIF677 678 c_u_m_l(k) = c_u_m_l(k) + c_u(k,i)679 c_v_m_l(k) = c_v_m_l(k) + c_v(k,i)680 c_w_m_l(k) = c_w_m_l(k) + c_w(k,i)681 682 ENDDO683 ENDDO684 685 #if defined( __parallel )686 IF ( collective_wait ) CALL MPI_BARRIER( comm1dx, ierr )687 CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, &688 MPI_SUM, comm1dx, ierr )689 IF ( collective_wait ) CALL MPI_BARRIER( comm1dx, ierr )690 CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, &691 MPI_SUM, comm1dx, ierr )692 IF ( collective_wait ) CALL MPI_BARRIER( comm1dx, ierr )693 CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, &694 MPI_SUM, comm1dx, ierr )695 #else696 c_u_m = c_u_m_l697 c_v_m = c_v_m_l698 c_w_m = c_w_m_l699 #endif700 701 c_u_m = c_u_m / (nx+1)702 c_v_m = c_v_m / (nx+1)703 c_w_m = c_w_m / (nx+1)704 705 !706 !-- Save old timelevels for the next timestep707 IF ( intermediate_timestep_count == 1 ) THEN708 630 u_m_n(:,:,:) = u(:,ny-1:ny,:) 709 631 v_m_n(:,:,:) = v(:,ny-1:ny,:) … … 731 653 u_p(nzb,ny+1,:) = 0.0_wp 732 654 v_p(nzb,ny+1,:) = 0.0_wp 733 ELSE 655 ELSE 734 656 u_p(nzb,ny+1,:) = u_p(nzb+1,ny+1,:) 735 657 v_p(nzb,ny+1,:) = v_p(nzb+1,ny+1,:) … … 757 679 u_p(:,:,0) = u(:,:,1) 758 680 v_p(:,:,-1) = v(:,:,0) 759 w_p(:,:,-1) = w(:,:,0) 681 w_p(:,:,-1) = w(:,:,0) 760 682 ELSEIF ( .NOT. use_cmax ) THEN 761 683 … … 870 792 u_p(nzb,:,0) = 0.0_wp 871 793 v_p(nzb,:,-1) = 0.0_wp 872 ELSE 794 ELSE 873 795 u_p(nzb,:,0) = u_p(nzb+1,:,0) 874 796 v_p(nzb,:,-1) = v_p(nzb+1,:,-1) … … 896 818 u_p(:,:,nx+1) = u(:,:,nx) 897 819 v_p(:,:,nx+1) = v(:,:,nx) 898 w_p(:,:,nx+1) = w(:,:,nx) 820 w_p(:,:,nx+1) = w(:,:,nx) 899 821 ELSEIF ( .NOT. use_cmax ) THEN 900 822 … … 1009 931 u_p(nzb,:,nx+1) = 0.0_wp 1010 932 v_p(nzb,:,nx+1) = 0.0_wp 1011 ELSE 933 ELSE 1012 934 u_p(nzb,:,nx+1) = u_p(nzb+1,:,nx+1) 1013 935 v_p(nzb,:,nx+1) = v_p(nzb+1,:,nx+1)
Note: See TracChangeset
for help on using the changeset viewer.
