Changeset 4102 for palm/trunk/SOURCE/boundary_conds.f90
- Timestamp:
- Jul 17, 2019 4:00:03 PM (5 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/boundary_conds.f90 (modified) (32 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/boundary_conds.f90
r4087 r4102 25 25 ! ----------------- 26 26 ! $Id$ 27 ! - Revise setting for boundary conditions at horizontal walls, use the offset 28 ! index that belongs to the data structure instead of pre-calculating 29 ! the offset index for each facing. 30 ! - Set boundary conditions for bulk-cloud quantities also at downward facing 31 ! surfaces 32 ! 33 ! 4087 2019-07-11 11:35:23Z gronemeier 27 34 ! Add comment 28 35 ! … … 224 231 USE arrays_3d, & 225 232 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, & 226 d iss, diss_p, dzu, e_p, nc_p, nr_p, pt, pt_init, pt_p, q,&233 dzu, nc_p, nr_p, pt, pt_init, pt_p, q, & 227 234 q_p, qc_p, qr_p, s, s_p, sa, sa_p, u, u_init, u_m_l, u_m_n, & 228 235 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, & 229 236 w, w_p, w_m_l, w_m_n, w_m_r, w_m_s 230 237 231 USE basic_constants_and_equations_mod, &232 ONLY: kappa233 234 238 USE bulk_cloud_model_mod, & 235 239 ONLY: bulk_cloud_model, microphysics_morrison, microphysics_seifert … … 239 243 240 244 USE control_parameters, & 241 ONLY: air_chemistry, bc_dirichlet_l, bc_dirichlet_n, bc_dirichlet_r,&245 ONLY: air_chemistry, bc_dirichlet_l, & 242 246 bc_dirichlet_s, bc_radiation_l, bc_radiation_n, bc_radiation_r, & 243 247 bc_radiation_s, bc_pt_t_val, bc_q_t_val, bc_s_t_val, & 244 child_domain, co nstant_diffusion, coupling_mode, dt_3d,&248 child_domain, coupling_mode, dt_3d, & 245 249 humidity, ibc_pt_b, ibc_pt_t, ibc_q_b, ibc_q_t, ibc_s_b, & 246 250 ibc_s_t, ibc_uv_b, ibc_uv_t, intermediate_timestep_count, & 247 251 nesting_offline, neutral, nudging, ocean_mode, passive_scalar, & 248 rans_mode, rans_tke_e,tsc, salsa, use_cmax252 tsc, salsa, use_cmax 249 253 250 254 USE grid_variables, & … … 262 266 263 267 USE pmc_interface, & 264 ONLY : nesting_mode , rans_mode_parent268 ONLY : nesting_mode 265 269 266 270 USE salsa_mod, & … … 268 272 269 273 USE surface_mod, & 270 ONLY : bc_h, surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, & 271 surf_usm_h, surf_usm_v 274 ONLY : bc_h 272 275 273 276 USE turbulence_closure_mod, & 274 ONLY: c_0277 ONLY: tcm_boundary_conds 275 278 276 279 IMPLICIT NONE … … 279 282 INTEGER(iwp) :: j !< grid index y direction 280 283 INTEGER(iwp) :: k !< grid index z direction 281 INTEGER(iwp) :: kb !< variable to set respective boundary value, depends on facing.282 284 INTEGER(iwp) :: l !< running index boundary type, for up- and downward-facing walls 283 285 INTEGER(iwp) :: m !< running index surface elements … … 297 299 !-- of downward-facing surfaces. 298 300 DO l = 0, 1 299 !300 !-- Set kb, for upward-facing surfaces value at topography top (k-1) is set,301 !-- for downward-facing surfaces at topography bottom (k+1).302 kb = MERGE( -1, 1, l == 0 )303 301 !$OMP PARALLEL DO PRIVATE( i, j, k ) 304 302 !$ACC PARALLEL LOOP PRIVATE(i, j, k) & … … 308 306 j = bc_h(l)%j(m) 309 307 k = bc_h(l)%k(m) 310 w_p(k+ kb,j,i) = 0.0_wp308 w_p(k+bc_h(l)%koff,j,i) = 0.0_wp 311 309 ENDDO 312 310 ENDDO … … 341 339 IF ( ibc_pt_b == 0 ) THEN 342 340 DO l = 0, 1 343 !344 !-- Set kb, for upward-facing surfaces value at topography top (k-1)345 !-- is set, for downward-facing surfaces at topography bottom (k+1).346 kb = MERGE( -1, 1, l == 0 )347 341 !$OMP PARALLEL DO PRIVATE( i, j, k ) 348 342 DO m = 1, bc_h(l)%ns … … 350 344 j = bc_h(l)%j(m) 351 345 k = bc_h(l)%k(m) 352 pt_p(k+ kb,j,i) = pt(k+kb,j,i)346 pt_p(k+bc_h(l)%koff,j,i) = pt(k+bc_h(l)%koff,j,i) 353 347 ENDDO 354 348 ENDDO … … 357 351 ELSEIF ( ibc_pt_b == 1 ) THEN 358 352 DO l = 0, 1 359 !360 !-- Set kb, for upward-facing surfaces value at topography top (k-1)361 !-- is set, for downward-facing surfaces at topography bottom (k+1).362 kb = MERGE( -1, 1, l == 0 )363 353 !$OMP PARALLEL DO PRIVATE( i, j, k ) 364 354 !$ACC PARALLEL LOOP PRIVATE(i, j, k) & … … 368 358 j = bc_h(l)%j(m) 369 359 k = bc_h(l)%k(m) 370 pt_p(k+ kb,j,i) = pt_p(k,j,i)360 pt_p(k+bc_h(l)%koff,j,i) = pt_p(k,j,i) 371 361 ENDDO 372 362 ENDDO … … 391 381 ENDIF 392 382 ENDIF 393 394 !395 !-- Boundary conditions for TKE.396 !-- Generally Neumann conditions with de/dz=0 are assumed.397 IF ( .NOT. constant_diffusion ) THEN398 399 IF ( .NOT. rans_mode ) THEN400 DO l = 0, 1401 !402 !-- Set kb, for upward-facing surfaces value at topography top (k-1) is set,403 !-- for downward-facing surfaces at topography bottom (k+1).404 kb = MERGE( -1, 1, l == 0 )405 !$OMP PARALLEL DO PRIVATE( i, j, k )406 !$ACC PARALLEL LOOP PRIVATE(i, j, k) &407 !$ACC PRESENT(bc_h, e_p)408 DO m = 1, bc_h(l)%ns409 i = bc_h(l)%i(m)410 j = bc_h(l)%j(m)411 k = bc_h(l)%k(m)412 e_p(k+kb,j,i) = e_p(k,j,i)413 ENDDO414 ENDDO415 ELSE416 !417 !-- Use wall function within constant-flux layer418 !-- Note, grid points listed in bc_h are not included in any calculations in RANS mode and419 !-- are therefore not set here.420 !421 !-- Upward-facing surfaces422 !-- Default surfaces423 DO m = 1, surf_def_h(0)%ns424 i = surf_def_h(0)%i(m)425 j = surf_def_h(0)%j(m)426 k = surf_def_h(0)%k(m)427 e_p(k,j,i) = ( surf_def_h(0)%us(m) / c_0 )**2428 ENDDO429 !430 !-- Natural surfaces431 DO m = 1, surf_lsm_h%ns432 i = surf_lsm_h%i(m)433 j = surf_lsm_h%j(m)434 k = surf_lsm_h%k(m)435 e_p(k,j,i) = ( surf_lsm_h%us(m) / c_0 )**2436 ENDDO437 !438 !-- Urban surfaces439 DO m = 1, surf_usm_h%ns440 i = surf_usm_h%i(m)441 j = surf_usm_h%j(m)442 k = surf_usm_h%k(m)443 e_p(k,j,i) = ( surf_usm_h%us(m) / c_0 )**2444 ENDDO445 !446 !-- Vertical surfaces447 DO l = 0, 3448 !449 !-- Default surfaces450 DO m = 1, surf_def_v(l)%ns451 i = surf_def_v(l)%i(m)452 j = surf_def_v(l)%j(m)453 k = surf_def_v(l)%k(m)454 e_p(k,j,i) = ( surf_def_v(l)%us(m) / c_0 )**2455 ENDDO456 !457 !-- Natural surfaces458 DO m = 1, surf_lsm_v(l)%ns459 i = surf_lsm_v(l)%i(m)460 j = surf_lsm_v(l)%j(m)461 k = surf_lsm_v(l)%k(m)462 e_p(k,j,i) = ( surf_lsm_v(l)%us(m) / c_0 )**2463 ENDDO464 !465 !-- Urban surfaces466 DO m = 1, surf_usm_v(l)%ns467 i = surf_usm_v(l)%i(m)468 j = surf_usm_v(l)%j(m)469 k = surf_usm_v(l)%k(m)470 e_p(k,j,i) = ( surf_usm_v(l)%us(m) / c_0 )**2471 ENDDO472 ENDDO473 ENDIF474 475 IF ( .NOT. child_domain ) THEN476 !$ACC KERNELS PRESENT(e_p)477 e_p(nzt+1,:,:) = e_p(nzt,:,:)478 !$ACC END KERNELS479 !480 !-- Nesting case: if parent operates in RANS mode and child in LES mode,481 !-- no TKE is transfered. This case, set Neumann conditions at lateral and482 !-- top child boundaries.483 !-- If not ( both either in RANS or in LES mode ), TKE boundary condition484 !-- is treated in the nesting.485 ELSE486 487 IF ( rans_mode_parent .AND. .NOT. rans_mode ) THEN488 489 e_p(nzt+1,:,:) = e_p(nzt,:,:)490 IF ( bc_dirichlet_l ) e_p(:,:,nxl-1) = e_p(:,:,nxl)491 IF ( bc_dirichlet_r ) e_p(:,:,nxr+1) = e_p(:,:,nxr)492 IF ( bc_dirichlet_s ) e_p(:,nys-1,:) = e_p(:,nys,:)493 IF ( bc_dirichlet_n ) e_p(:,nyn+1,:) = e_p(:,nyn,:)494 495 ENDIF496 ENDIF497 ENDIF498 499 !500 !-- Boundary conditions for TKE dissipation rate.501 IF ( rans_tke_e ) THEN502 !503 !-- Use wall function within constant-flux layer504 !-- Upward-facing surfaces505 !-- Default surfaces506 DO m = 1, surf_def_h(0)%ns507 i = surf_def_h(0)%i(m)508 j = surf_def_h(0)%j(m)509 k = surf_def_h(0)%k(m)510 diss_p(k,j,i) = surf_def_h(0)%us(m)**3 &511 / ( kappa * surf_def_h(0)%z_mo(m) )512 ENDDO513 !514 !-- Natural surfaces515 DO m = 1, surf_lsm_h%ns516 i = surf_lsm_h%i(m)517 j = surf_lsm_h%j(m)518 k = surf_lsm_h%k(m)519 diss_p(k,j,i) = surf_lsm_h%us(m)**3 &520 / ( kappa * surf_lsm_h%z_mo(m) )521 ENDDO522 !523 !-- Urban surfaces524 DO m = 1, surf_usm_h%ns525 i = surf_usm_h%i(m)526 j = surf_usm_h%j(m)527 k = surf_usm_h%k(m)528 diss_p(k,j,i) = surf_usm_h%us(m)**3 &529 / ( kappa * surf_usm_h%z_mo(m) )530 ENDDO531 !532 !-- Vertical surfaces533 DO l = 0, 3534 !535 !-- Default surfaces536 DO m = 1, surf_def_v(l)%ns537 i = surf_def_v(l)%i(m)538 j = surf_def_v(l)%j(m)539 k = surf_def_v(l)%k(m)540 diss_p(k,j,i) = surf_def_v(l)%us(m)**3 &541 / ( kappa * surf_def_v(l)%z_mo(m) )542 ENDDO543 !544 !-- Natural surfaces545 DO m = 1, surf_lsm_v(l)%ns546 i = surf_lsm_v(l)%i(m)547 j = surf_lsm_v(l)%j(m)548 k = surf_lsm_v(l)%k(m)549 diss_p(k,j,i) = surf_lsm_v(l)%us(m)**3 &550 / ( kappa * surf_lsm_v(l)%z_mo(m) )551 ENDDO552 !553 !-- Urban surfaces554 DO m = 1, surf_usm_v(l)%ns555 i = surf_usm_v(l)%i(m)556 j = surf_usm_v(l)%j(m)557 k = surf_usm_v(l)%k(m)558 diss_p(k,j,i) = surf_usm_v(l)%us(m)**3 &559 / ( kappa * surf_usm_v(l)%z_mo(m) )560 ENDDO561 ENDDO562 !563 !-- Limit change of diss to be between -90% and +100%. Also, set an absolute564 !-- minimum value565 DO i = nxl, nxr566 DO j = nys, nyn567 DO k = nzb, nzt+1568 diss_p(k,j,i) = MAX( MIN( diss_p(k,j,i), &569 2.0_wp * diss(k,j,i) ), &570 0.1_wp * diss(k,j,i), &571 0.0001_wp )572 ENDDO573 ENDDO574 ENDDO575 576 IF ( .NOT. child_domain ) THEN577 diss_p(nzt+1,:,:) = diss_p(nzt,:,:)578 ENDIF579 ENDIF580 581 383 ! 582 384 !-- Boundary conditions for salinity … … 586 388 !-- given. 587 389 DO l = 0, 1 588 !589 !-- Set kb, for upward-facing surfaces value at topography top (k-1) is set,590 !-- for downward-facing surfaces at topography bottom (k+1).591 kb = MERGE( -1, 1, l == 0 )592 390 !$OMP PARALLEL DO PRIVATE( i, j, k ) 593 391 DO m = 1, bc_h(l)%ns … … 595 393 j = bc_h(l)%j(m) 596 394 k = bc_h(l)%k(m) 597 sa_p(k+ kb,j,i) = sa_p(k,j,i)395 sa_p(k+bc_h(l)%koff,j,i) = sa_p(k,j,i) 598 396 ENDDO 599 397 ENDDO … … 620 418 621 419 DO l = 0, 1 622 !623 !-- Set kb, for upward-facing surfaces value at topography top (k-1) is set,624 !-- for downward-facing surfaces at topography bottom (k+1).625 kb = MERGE( -1, 1, l == 0 )626 420 !$OMP PARALLEL DO PRIVATE( i, j, k ) 627 421 DO m = 1, bc_h(l)%ns … … 629 423 j = bc_h(l)%j(m) 630 424 k = bc_h(l)%k(m) 631 q_p(k+ kb,j,i) = q(k+kb,j,i)425 q_p(k+bc_h(l)%koff,j,i) = q(k+bc_h(l)%koff,j,i) 632 426 ENDDO 633 427 ENDDO … … 636 430 637 431 DO l = 0, 1 638 !639 !-- Set kb, for upward-facing surfaces value at topography top (k-1) is set,640 !-- for downward-facing surfaces at topography bottom (k+1).641 kb = MERGE( -1, 1, l == 0 )642 432 !$OMP PARALLEL DO PRIVATE( i, j, k ) 643 433 DO m = 1, bc_h(l)%ns … … 645 435 j = bc_h(l)%j(m) 646 436 k = bc_h(l)%k(m) 647 q_p(k+ kb,j,i) = q_p(k,j,i)437 q_p(k+bc_h(l)%koff,j,i) = q_p(k,j,i) 648 438 ENDDO 649 439 ENDDO … … 662 452 !-- Run loop over all non-natural and natural walls. Note, in wall-datatype 663 453 !-- the k coordinate belongs to the atmospheric grid point, therefore, set 664 !-- qr_p and nr_p at k-1 454 !-- qr_p and nr_p at upward (k-1) and downward-facing (k+1) walls 455 DO l = 0, 1 665 456 !$OMP PARALLEL DO PRIVATE( i, j, k ) 666 DO m = 1, bc_h(0)%ns 667 i = bc_h(0)%i(m) 668 j = bc_h(0)%j(m) 669 k = bc_h(0)%k(m) 670 qc_p(k-1,j,i) = 0.0_wp 671 nc_p(k-1,j,i) = 0.0_wp 457 DO m = 1, bc_h(l)%ns 458 i = bc_h(l)%i(m) 459 j = bc_h(l)%j(m) 460 k = bc_h(l)%k(m) 461 qc_p(k+bc_h(l)%koff,j,i) = 0.0_wp 462 nc_p(k+bc_h(l)%koff,j,i) = 0.0_wp 463 ENDDO 672 464 ENDDO 673 465 ! … … 683 475 !-- Run loop over all non-natural and natural walls. Note, in wall-datatype 684 476 !-- the k coordinate belongs to the atmospheric grid point, therefore, set 685 !-- qr_p and nr_p at k-1 477 !-- qr_p and nr_p at upward (k-1) and downward-facing (k+1) walls 478 DO l = 0, 1 686 479 !$OMP PARALLEL DO PRIVATE( i, j, k ) 687 DO m = 1, bc_h(0)%ns 688 i = bc_h(0)%i(m) 689 j = bc_h(0)%j(m) 690 k = bc_h(0)%k(m) 691 qr_p(k-1,j,i) = 0.0_wp 692 nr_p(k-1,j,i) = 0.0_wp 480 DO m = 1, bc_h(l)%ns 481 i = bc_h(l)%i(m) 482 j = bc_h(l)%j(m) 483 k = bc_h(l)%k(m) 484 qr_p(k+bc_h(l)%koff,j,i) = 0.0_wp 485 nr_p(k+bc_h(l)%koff,j,i) = 0.0_wp 486 ENDDO 693 487 ENDDO 694 488 ! … … 711 505 712 506 DO l = 0, 1 713 !714 !-- Set kb, for upward-facing surfaces value at topography top (k-1) is set,715 !-- for downward-facing surfaces at topography bottom (k+1).716 kb = MERGE( -1, 1, l == 0 )717 507 !$OMP PARALLEL DO PRIVATE( i, j, k ) 718 508 DO m = 1, bc_h(l)%ns … … 720 510 j = bc_h(l)%j(m) 721 511 k = bc_h(l)%k(m) 722 s_p(k+ kb,j,i) = s(k+kb,j,i)512 s_p(k+bc_h(l)%koff,j,i) = s(k+bc_h(l)%koff,j,i) 723 513 ENDDO 724 514 ENDDO … … 727 517 728 518 DO l = 0, 1 729 !730 !-- Set kb, for upward-facing surfaces value at topography top (k-1) is set,731 !-- for downward-facing surfaces at topography bottom (k+1).732 kb = MERGE( -1, 1, l == 0 )733 519 !$OMP PARALLEL DO PRIVATE( i, j, k ) 734 520 DO m = 1, bc_h(l)%ns … … 736 522 j = bc_h(l)%j(m) 737 523 k = bc_h(l)%k(m) 738 s_p(k+ kb,j,i) = s_p(k,j,i)524 s_p(k+bc_h(l)%koff,j,i) = s_p(k,j,i) 739 525 ENDDO 740 526 ENDDO … … 750 536 ENDIF 751 537 752 ENDIF 538 ENDIF 539 ! 540 !-- Set boundary conditions for subgrid TKE and dissipation (RANS mode) 541 CALL tcm_boundary_conds 753 542 ! 754 543 !-- Top/bottom boundary conditions for chemical species … … 760 549 !-- version) these levels are handled as a prognostic level, boundary values 761 550 !-- have to be restored here. 762 !-- For the SGS-TKE, Neumann boundary conditions are used at the inflow.763 551 IF ( bc_dirichlet_s ) THEN 764 552 v_p(:,nys,:) = v_p(:,nys-1,:) 765 IF ( .NOT. constant_diffusion ) e_p(:,nys-1,:) = e_p(:,nys,:)766 ELSEIF ( bc_dirichlet_n ) THEN767 IF ( .NOT. constant_diffusion ) e_p(:,nyn+1,:) = e_p(:,nyn,:)768 553 ELSEIF ( bc_dirichlet_l ) THEN 769 554 u_p(:,:,nxl) = u_p(:,:,nxl-1) 770 IF ( .NOT. constant_diffusion ) e_p(:,:,nxl-1) = e_p(:,:,nxl)771 ELSEIF ( bc_dirichlet_r ) THEN772 IF ( .NOT. constant_diffusion ) e_p(:,:,nxr+1) = e_p(:,:,nxr)773 555 ENDIF 774 556 … … 777 559 !-- in case of nest boundaries. This must not be done in case of vertical nesting 778 560 !-- mode as in that case the lateral boundaries are actually cyclic. 561 !-- Lateral oundary conditions for TKE and dissipation are set 562 !-- in tcm_boundary_conds. 779 563 IF ( nesting_mode /= 'vertical' .OR. nesting_offline ) THEN 780 564 IF ( bc_dirichlet_s ) THEN … … 787 571 788 572 ! 789 !-- Lateral boundary conditions for scalar quantities at the outflow 573 !-- Lateral boundary conditions for scalar quantities at the outflow. 574 !-- Lateral oundary conditions for TKE and dissipation are set 575 !-- in tcm_boundary_conds. 790 576 IF ( bc_radiation_s ) THEN 791 577 pt_p(:,nys-1,:) = pt_p(:,nys,:) 792 IF ( .NOT. constant_diffusion ) e_p(:,nys-1,:) = e_p(:,nys,:)793 IF ( rans_tke_e ) diss_p(:,nys-1,:) = diss_p(:,nys,:)794 578 IF ( humidity ) THEN 795 579 q_p(:,nys-1,:) = q_p(:,nys,:) … … 806 590 ELSEIF ( bc_radiation_n ) THEN 807 591 pt_p(:,nyn+1,:) = pt_p(:,nyn,:) 808 IF ( .NOT. constant_diffusion ) e_p(:,nyn+1,:) = e_p(:,nyn,:)809 IF ( rans_tke_e ) diss_p(:,nyn+1,:) = diss_p(:,nyn,:)810 592 IF ( humidity ) THEN 811 593 q_p(:,nyn+1,:) = q_p(:,nyn,:) … … 822 604 ELSEIF ( bc_radiation_l ) THEN 823 605 pt_p(:,:,nxl-1) = pt_p(:,:,nxl) 824 IF ( .NOT. constant_diffusion ) e_p(:,:,nxl-1) = e_p(:,:,nxl)825 IF ( rans_tke_e ) diss_p(:,:,nxl-1) = diss_p(:,:,nxl)826 606 IF ( humidity ) THEN 827 607 q_p(:,:,nxl-1) = q_p(:,:,nxl) … … 838 618 ELSEIF ( bc_radiation_r ) THEN 839 619 pt_p(:,:,nxr+1) = pt_p(:,:,nxr) 840 IF ( .NOT. constant_diffusion ) e_p(:,:,nxr+1) = e_p(:,:,nxr)841 IF ( rans_tke_e ) diss_p(:,:,nxr+1) = diss_p(:,:,nxr)842 620 IF ( humidity ) THEN 843 621 q_p(:,:,nxr+1) = q_p(:,:,nxr)
Note: See TracChangeset
for help on using the changeset viewer.
