Changeset 4102 for palm/trunk/SOURCE
- Timestamp:
- Jul 17, 2019 4:00:03 PM (5 years ago)
- Location:
- palm/trunk/SOURCE
- Files:
-
- 7 edited
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/Makefile
r4070 r4102 25 25 # ----------------- 26 26 # $Id$ 27 # Remove dependency on pmc_interface for boundary_conds 28 # 29 # 4070 2019-07-03 13:51:40Z gronemeier 27 30 # Add new data output modules 28 31 # … … 802 805 radiation_model_mod.o 803 806 boundary_conds.o: \ 804 basic_constants_and_equations_mod.o \805 807 bulk_cloud_model_mod.o \ 806 808 chemistry_model_mod.o \ 807 809 mod_kinds.o \ 808 810 modules.o \ 809 pmc_interface_mod.o \810 811 salsa_mod.o \ 811 812 surface_mod.o \ -
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) -
palm/trunk/SOURCE/chemistry_model_mod.f90
r4080 r4102 27 27 ! ----------------- 28 28 ! $Id$ 29 ! Slightly revise setting of boundary conditions at horizontal walls, use 30 ! data-structure offset index instead of pre-calculate it for each facing 31 ! 32 ! 4080 2019-07-09 18:17:37Z suehring 29 33 ! Restore accidantly removed limitation to positive values 30 34 ! … … 804 808 INTEGER(iwp) :: j !< grid index y direction. 805 809 INTEGER(iwp) :: k !< grid index z direction. 806 INTEGER(iwp) :: kb !< variable to set respective boundary value, depends on facing.807 810 INTEGER(iwp) :: l !< running index boundary type, for up- and downward-facing walls. 808 811 INTEGER(iwp) :: m !< running index surface elements. … … 839 842 CASE ( 'set_bc_bottomtop' ) 840 843 ! 841 !-- Bottom boundary condtions for chemical species844 !-- Boundary condtions for chemical species at horizontal walls 842 845 DO lsp = 1, nspec 843 IF ( ibc_cs_b == 0 ) THEN 846 IF ( ibc_cs_b == 0 ) THEN 844 847 DO l = 0, 1 845 !846 !-- Set index kb: For upward-facing surfaces (l=0), kb=-1, i.e.847 !-- the chem_species(nspec)%conc_p value at the topography top (k-1)848 !-- is set; for downward-facing surfaces (l=1), kb=1, i.e. the849 !-- value at the topography bottom (k+1) is set.850 851 kb = MERGE( -1, 1, l == 0 )852 848 !$OMP PARALLEL DO PRIVATE( i, j, k ) 853 849 DO m = 1, bc_h(l)%ns … … 855 851 j = bc_h(l)%j(m) 856 852 k = bc_h(l)%k(m) 857 chem_species(lsp)%conc_p(k+kb,j,i) = chem_species(lsp)%conc(k+kb,j,i) 853 chem_species(lsp)%conc_p(k+bc_h(l)%koff,j,i) = & 854 chem_species(lsp)%conc(k+bc_h(l)%koff,j,i) 858 855 ENDDO 859 856 ENDDO … … 863 860 !-- in boundary_conds there is som extra loop over m here for passive tracer 864 861 DO l = 0, 1 865 kb = MERGE( -1, 1, l == 0 )866 862 !$OMP PARALLEL DO PRIVATE( i, j, k ) 867 863 DO m = 1, bc_h(l)%ns … … 869 865 j = bc_h(l)%j(m) 870 866 k = bc_h(l)%k(m) 871 chem_species(lsp)%conc_p(k+kb,j,i) = chem_species(lsp)%conc_p(k,j,i) 867 chem_species(lsp)%conc_p(k+bc_h(l)%koff,j,i) = & 868 chem_species(lsp)%conc_p(k,j,i) 872 869 873 870 ENDDO -
palm/trunk/SOURCE/salsa_mod.f90
r4079 r4102 26 26 ! ----------------- 27 27 ! $Id$ 28 ! Slightly revise setting of boundary conditions at horizontal walls, use 29 ! data-structure offset index instead of pre-calculate it for each facing 30 ! 31 ! 4079 2019-07-09 18:04:41Z suehring 28 32 ! Application of monotonic flux limiter for the vertical scalar advection 29 33 ! up to the topography top (only for the cache-optimized version at the … … 8017 8021 INTEGER(iwp) :: j !< grid index y direction 8018 8022 INTEGER(iwp) :: k !< grid index y direction 8019 INTEGER(iwp) :: kb !< variable to set respective boundary value, depends on facing.8020 8023 INTEGER(iwp) :: l !< running index boundary type, for up- and downward-facing walls 8021 8024 INTEGER(iwp) :: m !< running index surface elements … … 8030 8033 !-- belongs to the atmospheric grid point, therefore, set s_p at k-1 8031 8034 DO l = 0, 1 8032 !8033 !-- Set kb, for upward-facing surfaces value at topography top (k-1) is8034 !-- set, for downward-facing surfaces at topography bottom (k+1)8035 kb = MERGE ( -1, 1, l == 0 )8036 8035 !$OMP PARALLEL PRIVATE( ib, ic, icc, ig, i, j, k ) 8037 8036 !$OMP DO … … 8043 8042 8044 8043 DO ib = 1, nbins_aerosol 8045 aerosol_number(ib)%conc_p(k+kb,j,i) = aerosol_number(ib)%conc(k+kb,j,i) 8044 aerosol_number(ib)%conc_p(k+bc_h(l)%koff,j,i) = & 8045 aerosol_number(ib)%conc(k+bc_h(l)%koff,j,i) 8046 8046 DO ic = 1, ncomponents_mass 8047 8047 icc = ( ic - 1 ) * nbins_aerosol + ib 8048 aerosol_mass(icc)%conc_p(k+kb,j,i) = aerosol_mass(icc)%conc(k+kb,j,i) 8048 aerosol_mass(icc)%conc_p(k+bc_h(l)%koff,j,i) = & 8049 aerosol_mass(icc)%conc(k+bc_h(l)%koff,j,i) 8049 8050 ENDDO 8050 8051 ENDDO 8051 8052 IF ( .NOT. salsa_gases_from_chem ) THEN 8052 8053 DO ig = 1, ngases_salsa 8053 salsa_gas(ig)%conc_p(k+kb,j,i) = salsa_gas(ig)%conc(k+kb,j,i) 8054 salsa_gas(ig)%conc_p(k+bc_h(l)%koff,j,i) = & 8055 salsa_gas(ig)%conc(k+bc_h(l)%koff,j,i) 8054 8056 ENDDO 8055 8057 ENDIF … … 8063 8065 8064 8066 DO l = 0, 1 8065 !8066 !-- Set kb, for upward-facing surfaces value at topography top (k-1) is8067 !-- set, for downward-facing surfaces at topography bottom (k+1)8068 kb = MERGE( -1, 1, l == 0 )8069 8067 !$OMP PARALLEL PRIVATE( ib, ic, icc, ig, i, j, k ) 8070 8068 !$OMP DO … … 8076 8074 8077 8075 DO ib = 1, nbins_aerosol 8078 aerosol_number(ib)%conc_p(k+kb,j,i) = aerosol_number(ib)%conc_p(k,j,i) 8076 aerosol_number(ib)%conc_p(k+bc_h(l)%koff,j,i) = & 8077 aerosol_number(ib)%conc_p(k,j,i) 8079 8078 DO ic = 1, ncomponents_mass 8080 8079 icc = ( ic - 1 ) * nbins_aerosol + ib 8081 aerosol_mass(icc)%conc_p(k+kb,j,i) = aerosol_mass(icc)%conc_p(k,j,i) 8080 aerosol_mass(icc)%conc_p(k+bc_h(l)%koff,j,i) = & 8081 aerosol_mass(icc)%conc_p(k,j,i) 8082 8082 ENDDO 8083 8083 ENDDO 8084 8084 IF ( .NOT. salsa_gases_from_chem ) THEN 8085 8085 DO ig = 1, ngases_salsa 8086 salsa_gas(ig)%conc_p(k+kb,j,i) = salsa_gas(ig)%conc_p(k,j,i) 8086 salsa_gas(ig)%conc_p(k+bc_h(l)%koff,j,i) = & 8087 salsa_gas(ig)%conc_p(k,j,i) 8087 8088 ENDDO 8088 8089 ENDIF -
palm/trunk/SOURCE/surface_mod.f90
r3943 r4102 26 26 ! ----------------- 27 27 ! $Id$ 28 ! - Revise initialization of the boundary data structure 29 ! - Add new data structure to set boundary conditions at vertical walls 30 ! 31 ! 3943 2019-05-02 09:50:41Z maronga 28 32 ! Removed qsws_eb as it is no longer needed. 29 33 ! … … 306 310 !-- are applied 307 311 TYPE bc_type 308 309 INTEGER(iwp) :: ns !< number of surface elements on the PE 310 311 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: i !< x-index linking to the PALM 3D-grid 312 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: j !< y-index linking to the PALM 3D-grid 313 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: k !< z-index linking to the PALM 3D-grid 312 INTEGER(iwp) :: ioff !< offset value in x-direction, used to determine index of surface element 313 INTEGER(iwp) :: joff !< offset value in y-direction, used to determine index of surface element 314 INTEGER(iwp) :: koff !< offset value in z-direction, used to determine index of surface element 315 INTEGER(iwp) :: ns !< number of surface elements on the PE 316 317 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: i !< x-index linking to the PALM 3D-grid 318 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: j !< y-index linking to the PALM 3D-grid 319 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: k !< z-index linking to the PALM 3D-grid 314 320 315 321 INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: start_index !< start index within surface data type for given (j,i) … … 600 606 601 607 TYPE (bc_type), DIMENSION(0:1) :: bc_h !< boundary condition data type, horizontal upward- and downward facing surfaces 608 TYPE (bc_type), DIMENSION(0:3) :: bc_v !< boundary condition data type, vertical surfaces 602 609 603 610 TYPE (surf_type), DIMENSION(0:2), TARGET :: surf_def_h !< horizontal default surfaces (Up, Down, and Top) … … 668 675 ! 669 676 !-- Public variables 670 PUBLIC bc_h, ind_pav_green, ind_veg_wall, ind_wat_win, ns_h_on_file, ns_v_on_file, surf_def_h,&671 surf_def_ v, surf_lsm_h, surf_lsm_v, surf_usm_h, surf_usm_v, surf_type,&677 PUBLIC bc_h, bc_v, ind_pav_green, ind_veg_wall, ind_wat_win, ns_h_on_file, ns_v_on_file, & 678 surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, surf_usm_v, surf_type, & 672 679 vertical_surfaces_exist, surf_bulk_cloud_model, surf_microphysics_morrison, & 673 680 surf_microphysics_seifert … … 687 694 ! Description: 688 695 ! ------------ 689 !> Initialize data type for setting boundary conditions at horizontal surfaces. 696 !> Initialize data type for setting boundary conditions at horizontal and 697 !> vertical surfaces. 690 698 !------------------------------------------------------------------------------! 691 699 SUBROUTINE init_bc … … 696 704 INTEGER(iwp) :: j !< loop index along y-direction 697 705 INTEGER(iwp) :: k !< loop index along y-direction 698 699 INTEGER(iwp), DIMENSION(0:1) :: num_h !< number of surfaces 700 INTEGER(iwp), DIMENSION(0:1) :: num_h_kji !< number of surfaces for (j,i)-grid point 701 INTEGER(iwp), DIMENSION(0:1) :: start_index_h !< local start index 702 703 ! 704 !-- First of all, count the number of upward- and downward-facing surfaces 705 num_h = 0 706 DO i = nxlg, nxrg 707 DO j = nysg, nyng 708 DO k = nzb+1, nzt 709 ! 710 !-- Check if current gridpoint belongs to the atmosphere 711 IF ( BTEST( wall_flags_0(k,j,i), 0 ) ) THEN 712 ! 713 !-- Upward-facing 714 IF ( .NOT. BTEST( wall_flags_0(k-1,j,i), 0 ) ) & 715 num_h(0) = num_h(0) + 1 716 ! 717 !-- Downward-facing 718 IF ( .NOT. BTEST( wall_flags_0(k+1,j,i), 0 ) ) & 719 num_h(1) = num_h(1) + 1 720 ENDIF 706 INTEGER(iwp) :: l !< running index for differently aligned surfaces 707 708 INTEGER(iwp), DIMENSION(0:1) :: num_h !< number of horizontal surfaces on subdomain 709 INTEGER(iwp), DIMENSION(0:1) :: num_h_kji !< number of horizontal surfaces at (j,i)-grid point 710 INTEGER(iwp), DIMENSION(0:1) :: start_index_h !< local start index of horizontal surface elements 711 712 INTEGER(iwp), DIMENSION(0:3) :: num_v !< number of vertical surfaces on subdomain 713 INTEGER(iwp), DIMENSION(0:3) :: num_v_kji !< number of vertical surfaces at (j,i)-grid point 714 INTEGER(iwp), DIMENSION(0:3) :: start_index_v !< local start index of vertical surface elements 715 ! 716 !-- Set offset indices, i.e. index difference between surface element and 717 !-- surface-bounded grid point. 718 !-- Horizontal surfaces - no horizontal offsets 719 bc_h(:)%ioff = 0 720 bc_h(:)%joff = 0 721 ! 722 !-- Horizontal surfaces, upward facing (0) and downward facing (1) 723 bc_h(0)%koff = -1 724 bc_h(1)%koff = 1 725 ! 726 !-- Vertical surfaces - no vertical offset 727 bc_v(0:3)%koff = 0 728 ! 729 !-- North- and southward facing - no offset in x 730 bc_v(0:1)%ioff = 0 731 ! 732 !-- Northward facing offset in y 733 bc_v(0)%joff = -1 734 ! 735 !-- Southward facing offset in y 736 bc_v(1)%joff = 1 737 ! 738 !-- East- and westward facing - no offset in y 739 bc_v(2:3)%joff = 0 740 ! 741 !-- Eastward facing offset in x 742 bc_v(2)%ioff = -1 743 ! 744 !-- Westward facing offset in y 745 bc_v(3)%ioff = 1 746 ! 747 !-- Initialize data structure for horizontal surfaces, i.e. count the number 748 !-- of surface elements, allocate and initialize the respective index arrays, 749 !-- and set the respective start and end indices at each (j,i)-location. 750 DO l = 0, 1 751 ! 752 !-- Count the number of upward- and downward-facing surfaces on subdomain 753 num_h(l) = 0 754 DO i = nxlg, nxrg 755 DO j = nysg, nyng 756 DO k = nzb+1, nzt 757 ! 758 !-- Check if current gridpoint belongs to the atmosphere 759 IF ( BTEST( wall_flags_0(k,j,i), 0 ) ) THEN 760 IF ( .NOT. BTEST( wall_flags_0(k+bc_h(l)%koff, & 761 j+bc_h(l)%joff, & 762 i+bc_h(l)%ioff), 0 ) ) & 763 num_h(l) = num_h(l) + 1 764 ENDIF 765 ENDDO 766 ENDDO 767 ENDDO 768 ! 769 !-- Save the number of horizontal surface elements 770 bc_h(l)%ns = num_h(l) 771 ! 772 !-- ALLOCATE arrays for horizontal surfaces 773 ALLOCATE( bc_h(l)%i(1:bc_h(l)%ns) ) 774 ALLOCATE( bc_h(l)%j(1:bc_h(l)%ns) ) 775 ALLOCATE( bc_h(l)%k(1:bc_h(l)%ns) ) 776 ALLOCATE( bc_h(l)%start_index(nysg:nyng,nxlg:nxrg) ) 777 ALLOCATE( bc_h(l)%end_index(nysg:nyng,nxlg:nxrg) ) 778 bc_h(l)%start_index = 1 779 bc_h(l)%end_index = 0 780 781 num_h(l) = 1 782 start_index_h(l) = 1 783 DO i = nxlg, nxrg 784 DO j = nysg, nyng 785 786 num_h_kji(l) = 0 787 DO k = nzb+1, nzt 788 ! 789 !-- Check if current gridpoint belongs to the atmosphere 790 IF ( BTEST( wall_flags_0(k,j,i), 0 ) ) THEN 791 ! 792 !-- Upward-facing 793 IF ( .NOT. BTEST( wall_flags_0(k+bc_h(l)%koff, & 794 j+bc_h(l)%joff, & 795 i+bc_h(l)%ioff), 0 ) & 796 ) THEN 797 bc_h(l)%i(num_h(l)) = i 798 bc_h(l)%j(num_h(l)) = j 799 bc_h(l)%k(num_h(l)) = k 800 num_h_kji(l) = num_h_kji(l) + 1 801 num_h(l) = num_h(l) + 1 802 ENDIF 803 ENDIF 804 ENDDO 805 bc_h(l)%start_index(j,i) = start_index_h(l) 806 bc_h(l)%end_index(j,i) = bc_h(l)%start_index(j,i) + & 807 num_h_kji(l) - 1 808 start_index_h(l) = bc_h(l)%end_index(j,i) + 1 721 809 ENDDO 722 810 ENDDO 723 811 ENDDO 724 ! 725 !-- Save the number of surface elements 726 bc_h(0)%ns = num_h(0) 727 bc_h(1)%ns = num_h(1) 728 ! 729 !-- ALLOCATE data type variables 730 !-- Upward facing 731 ALLOCATE( bc_h(0)%i(1:bc_h(0)%ns) ) 732 ALLOCATE( bc_h(0)%j(1:bc_h(0)%ns) ) 733 ALLOCATE( bc_h(0)%k(1:bc_h(0)%ns) ) 734 ALLOCATE( bc_h(0)%start_index(nysg:nyng,nxlg:nxrg) ) 735 ALLOCATE( bc_h(0)%end_index(nysg:nyng,nxlg:nxrg) ) 736 bc_h(0)%start_index = 1 737 bc_h(0)%end_index = 0 738 ! 739 !-- Downward facing 740 ALLOCATE( bc_h(1)%i(1:bc_h(1)%ns) ) 741 ALLOCATE( bc_h(1)%j(1:bc_h(1)%ns) ) 742 ALLOCATE( bc_h(1)%k(1:bc_h(1)%ns) ) 743 ALLOCATE( bc_h(1)%start_index(nysg:nyng,nxlg:nxrg) ) 744 ALLOCATE( bc_h(1)%end_index(nysg:nyng,nxlg:nxrg) ) 745 bc_h(1)%start_index = 1 746 bc_h(1)%end_index = 0 747 ! 748 !-- Store the respective indices on data type 749 num_h(0:1) = 1 750 start_index_h(0:1) = 1 751 DO i = nxlg, nxrg 752 DO j = nysg, nyng 753 754 num_h_kji(0:1) = 0 755 DO k = nzb+1, nzt 756 ! 757 !-- Check if current gridpoint belongs to the atmosphere 758 IF ( BTEST( wall_flags_0(k,j,i), 0 ) ) THEN 759 ! 760 !-- Upward-facing 761 IF ( .NOT. BTEST( wall_flags_0(k-1,j,i), 0 ) ) THEN 762 bc_h(0)%i(num_h(0)) = i 763 bc_h(0)%j(num_h(0)) = j 764 bc_h(0)%k(num_h(0)) = k 765 num_h_kji(0) = num_h_kji(0) + 1 766 num_h(0) = num_h(0) + 1 812 813 ! 814 !-- Initialize data structure for vertical surfaces, i.e. count the number 815 !-- of surface elements, allocate and initialize the respective index arrays, 816 !-- and set the respective start and end indices at each (j,i)-location. 817 DO l = 0, 3 818 ! 819 !-- Count the number of upward- and downward-facing surfaces on subdomain 820 num_v(l) = 0 821 DO i = nxlg, nxrg 822 DO j = nysg, nyng 823 DO k = nzb+1, nzt 824 ! 825 !-- Check if current gridpoint belongs to the atmosphere 826 IF ( BTEST( wall_flags_0(k,j,i), 0 ) ) THEN 827 IF ( .NOT. BTEST( wall_flags_0(k+bc_v(l)%koff, & 828 j+bc_v(l)%joff, & 829 i+bc_v(l)%ioff), 0 ) ) & 830 num_v(l) = num_v(l) + 1 767 831 ENDIF 768 ! 769 !-- Downward-facing 770 IF ( .NOT. BTEST( wall_flags_0(k+1,j,i), 0 ) ) THEN 771 bc_h(1)%i(num_h(1)) = i 772 bc_h(1)%j(num_h(1)) = j 773 bc_h(1)%k(num_h(1)) = k 774 num_h_kji(1) = num_h_kji(1) + 1 775 num_h(1) = num_h(1) + 1 832 ENDDO 833 ENDDO 834 ENDDO 835 ! 836 !-- Save the number of horizontal surface elements 837 bc_v(l)%ns = num_v(l) 838 ! 839 !-- ALLOCATE arrays for horizontal surfaces 840 ALLOCATE( bc_v(l)%i(1:bc_v(l)%ns) ) 841 ALLOCATE( bc_v(l)%j(1:bc_v(l)%ns) ) 842 ALLOCATE( bc_v(l)%k(1:bc_v(l)%ns) ) 843 ALLOCATE( bc_v(l)%start_index(nysg:nyng,nxlg:nxrg) ) 844 ALLOCATE( bc_v(l)%end_index(nysg:nyng,nxlg:nxrg) ) 845 bc_v(l)%start_index = 1 846 bc_v(l)%end_index = 0 847 848 num_v(l) = 1 849 start_index_v(l) = 1 850 DO i = nxlg, nxrg 851 DO j = nysg, nyng 852 853 num_v_kji(l) = 0 854 DO k = nzb+1, nzt 855 ! 856 !-- Check if current gridpoint belongs to the atmosphere 857 IF ( BTEST( wall_flags_0(k,j,i), 0 ) ) THEN 858 ! 859 !-- Upward-facing 860 IF ( .NOT. BTEST( wall_flags_0(k+bc_v(l)%koff, & 861 j+bc_v(l)%joff, & 862 i+bc_v(l)%ioff), 0 ) & 863 ) THEN 864 bc_v(l)%i(num_v(l)) = i 865 bc_v(l)%j(num_v(l)) = j 866 bc_v(l)%k(num_v(l)) = k 867 num_v_kji(l) = num_v_kji(l) + 1 868 num_v(l) = num_v(l) + 1 869 ENDIF 776 870 ENDIF 777 ENDIF 871 ENDDO 872 bc_v(l)%start_index(j,i) = start_index_v(l) 873 bc_v(l)%end_index(j,i) = bc_v(l)%start_index(j,i) + & 874 num_v_kji(l) - 1 875 start_index_v(l) = bc_v(l)%end_index(j,i) + 1 778 876 ENDDO 779 bc_h(0)%start_index(j,i) = start_index_h(0)780 bc_h(0)%end_index(j,i) = bc_h(0)%start_index(j,i) + num_h_kji(0) - 1781 start_index_h(0) = bc_h(0)%end_index(j,i) + 1782 783 bc_h(1)%start_index(j,i) = start_index_h(1)784 bc_h(1)%end_index(j,i) = bc_h(1)%start_index(j,i) + num_h_kji(1) - 1785 start_index_h(1) = bc_h(1)%end_index(j,i) + 1786 877 ENDDO 787 878 ENDDO -
palm/trunk/SOURCE/turbulence_closure_mod.f90
r4048 r4102 25 25 ! ----------------- 26 26 ! $Id$ 27 ! - Modularize setting of boundary conditions for TKE and dissipation 28 ! - Neumann boundary condition for TKE at model top is set also in child domain 29 ! - Revise setting of Neumann boundary conditions at non-cyclic lateral 30 ! boundaries 31 ! - Bugfix, set Neumann boundary condition for TKE at vertical wall instead of 32 ! an implicit Dirichlet boundary condition which implied a sink of TKE 33 ! at vertical walls 34 ! 35 ! 4048 2019-06-21 21:00:21Z knoop 27 36 ! write out preprocessor directives; remove tailing whitespaces 28 37 ! … … 192 201 !> @todo Check for random disturbances 193 202 !> @note <Enter notes on the module> 194 !----------------------------------------------------------------------------- -!203 !-----------------------------------------------------------------------------! 195 204 MODULE turbulence_closure_mod 196 205 197 206 198 USE arrays_3d, 199 ONLY: diss, diss_1, diss_2, diss_3, diss_p, dzu, e, e_1, e_2, e_3, 200 e_p, kh, km, mean_inflow_profiles, prho, pt, tdiss_m, 207 USE arrays_3d, & 208 ONLY: diss, diss_1, diss_2, diss_3, diss_p, dzu, e, e_1, e_2, e_3, & 209 e_p, kh, km, mean_inflow_profiles, prho, pt, tdiss_m, & 201 210 te_m, tend, u, v, vpt, w 202 211 203 USE basic_constants_and_equations_mod, 212 USE basic_constants_and_equations_mod, & 204 213 ONLY: g, kappa, lv_d_cp, lv_d_rd, rd_d_rv 205 214 206 USE control_parameters, & 207 ONLY: constant_diffusion, dt_3d, e_init, humidity, & 208 initializing_actions, intermediate_timestep_count, & 209 intermediate_timestep_count_max, km_constant, & 210 les_dynamic, les_mw, ocean_mode, plant_canopy, prandtl_number, & 211 pt_reference, rans_mode, rans_tke_e, rans_tke_l, & 212 timestep_scheme, turbulence_closure, & 213 turbulent_inflow, use_upstream_for_tke, vpt_reference, & 215 USE control_parameters, & 216 ONLY: bc_dirichlet_l, & 217 bc_dirichlet_n, & 218 bc_dirichlet_r, & 219 bc_dirichlet_s, & 220 bc_radiation_l, & 221 bc_radiation_n, & 222 bc_radiation_r, & 223 bc_radiation_s, & 224 child_domain, & 225 constant_diffusion, dt_3d, e_init, humidity, & 226 initializing_actions, intermediate_timestep_count, & 227 intermediate_timestep_count_max, km_constant, & 228 les_dynamic, les_mw, ocean_mode, plant_canopy, prandtl_number, & 229 pt_reference, rans_mode, rans_tke_e, rans_tke_l, & 230 timestep_scheme, turbulence_closure, & 231 turbulent_inflow, use_upstream_for_tke, vpt_reference, & 214 232 ws_scheme_sca, current_timestep_number 215 233 216 USE advec_ws, 234 USE advec_ws, & 217 235 ONLY: advec_s_ws 218 236 219 USE advec_s_bc_mod, 237 USE advec_s_bc_mod, & 220 238 ONLY: advec_s_bc 221 239 222 USE advec_s_pw_mod, 240 USE advec_s_pw_mod, & 223 241 ONLY: advec_s_pw 224 242 225 USE advec_s_up_mod, 243 USE advec_s_up_mod, & 226 244 ONLY: advec_s_up 227 245 228 USE cpulog, 246 USE cpulog, & 229 247 ONLY: cpu_log, log_point_s 230 248 231 USE indices, 232 ONLY: nbgp, nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt, 249 USE indices, & 250 ONLY: nbgp, nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt, & 233 251 wall_flags_0 234 252 235 253 USE kinds 236 254 237 USE ocean_mod, 255 USE ocean_mod, & 238 256 ONLY: prho_reference 239 257 240 258 USE pegrid 241 259 242 USE plant_canopy_model_mod, 260 USE plant_canopy_model_mod, & 243 261 ONLY: pcm_tendency 244 262 245 USE statistics, 263 USE statistics, & 246 264 ONLY: hom, hom_sum, statistic_regions 247 265 266 USE surface_mod, & 267 ONLY: bc_h, & 268 bc_v, & 269 get_topography_top_index_ji, & 270 surf_def_h, & 271 surf_def_v, & 272 surf_lsm_h, & 273 surf_lsm_v, & 274 surf_usm_h, & 275 surf_usm_v 248 276 249 277 IMPLICIT NONE … … 271 299 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: l_wall !< near-wall mixing length 272 300 273 301 ! 302 !-- Public variables 274 303 PUBLIC c_0, rans_const_c, rans_const_sigma 275 304 … … 277 306 278 307 PRIVATE 279 280 281 PUBLIC & 282 tcm_check_parameters, & 283 tcm_check_data_output, & 284 tcm_define_netcdf_grid, & 285 tcm_init_arrays, & 286 tcm_init, & 287 tcm_actions, & 288 tcm_prognostic_equations, & 289 tcm_swap_timelevel, & 290 tcm_3d_data_averaging, & 291 tcm_data_output_2d, & 292 tcm_data_output_3d, & 308 ! 309 !-- Public subroutines 310 PUBLIC & 311 tcm_boundary_conds, & 312 tcm_check_parameters, & 313 tcm_check_data_output, & 314 tcm_define_netcdf_grid, & 315 tcm_init_arrays, & 316 tcm_init, & 317 tcm_actions, & 318 tcm_prognostic_equations, & 319 tcm_swap_timelevel, & 320 tcm_3d_data_averaging, & 321 tcm_data_output_2d, & 322 tcm_data_output_3d, & 293 323 tcm_diffusivities 294 324 295 325 ! 296 326 !-- PALM interfaces: 327 !-- Boundary conditions for subgrid TKE and dissipation 328 INTERFACE tcm_boundary_conds 329 MODULE PROCEDURE tcm_boundary_conds 330 END INTERFACE tcm_boundary_conds 331 ! 297 332 !-- Input parameter checks to be done in check_parameters 298 333 INTERFACE tcm_check_parameters … … 371 406 CONTAINS 372 407 408 !------------------------------------------------------------------------------! 409 ! Description: 410 ! ------------ 411 !> Check parameters routine for turbulence closure module. 412 !------------------------------------------------------------------------------! 413 SUBROUTINE tcm_boundary_conds 414 415 USE pmc_interface, & 416 ONLY : rans_mode_parent 417 418 IMPLICIT NONE 419 420 INTEGER(iwp) :: i !< grid index x direction 421 INTEGER(iwp) :: j !< grid index y direction 422 INTEGER(iwp) :: k !< grid index z direction 423 INTEGER(iwp) :: l !< running index boundary type, for up- and downward-facing walls 424 INTEGER(iwp) :: m !< running index surface elements 425 ! 426 !-- Boundary conditions for TKE. 427 IF ( .NOT. constant_diffusion ) THEN 428 ! 429 !-- In LES mode, Neumann conditions with de/x_i=0 are assumed at solid walls. 430 !-- Note, only TKE is prognostic in this case and dissipation is only 431 !-- a diagnostic quantity. 432 IF ( .NOT. rans_mode ) THEN 433 ! 434 !-- Horizontal walls, upward- and downward-facing 435 DO l = 0, 1 436 !$OMP PARALLEL DO PRIVATE( i, j, k ) 437 !$ACC PARALLEL LOOP PRIVATE(i, j, k) & 438 !$ACC PRESENT(bc_h, e_p) 439 DO m = 1, bc_h(l)%ns 440 i = bc_h(l)%i(m) 441 j = bc_h(l)%j(m) 442 k = bc_h(l)%k(m) 443 e_p(k+bc_h(l)%koff,j,i) = e_p(k,j,i) 444 ENDDO 445 ENDDO 446 ! 447 !-- Vertical walls 448 DO l = 0, 3 449 !$OMP PARALLEL DO PRIVATE( i, j, k ) 450 !$ACC PARALLEL LOOP PRIVATE(i, j, k) & 451 !$ACC PRESENT(bc_v, e_p) 452 DO m = 1, bc_v(l)%ns 453 i = bc_v(l)%i(m) 454 j = bc_v(l)%j(m) 455 k = bc_v(l)%k(m) 456 e_p(k,j+bc_v(l)%joff,i+bc_v(l)%ioff) = e_p(k,j,i) 457 ENDDO 458 ENDDO 459 ! 460 !-- In RANS mode, wall function is used as boundary condition for TKE 461 ELSE 462 ! 463 !-- Use wall function within constant-flux layer 464 !-- Note, grid points listed in bc_h are not included in any calculations in RANS mode and 465 !-- are therefore not set here. 466 ! 467 !-- Upward-facing surfaces 468 !-- Default surfaces 469 DO m = 1, surf_def_h(0)%ns 470 i = surf_def_h(0)%i(m) 471 j = surf_def_h(0)%j(m) 472 k = surf_def_h(0)%k(m) 473 e_p(k,j,i) = ( surf_def_h(0)%us(m) / c_0 )**2 474 ENDDO 475 ! 476 !-- Natural surfaces 477 DO m = 1, surf_lsm_h%ns 478 i = surf_lsm_h%i(m) 479 j = surf_lsm_h%j(m) 480 k = surf_lsm_h%k(m) 481 e_p(k,j,i) = ( surf_lsm_h%us(m) / c_0 )**2 482 ENDDO 483 ! 484 !-- Urban surfaces 485 DO m = 1, surf_usm_h%ns 486 i = surf_usm_h%i(m) 487 j = surf_usm_h%j(m) 488 k = surf_usm_h%k(m) 489 e_p(k,j,i) = ( surf_usm_h%us(m) / c_0 )**2 490 ENDDO 491 ! 492 !-- Vertical surfaces 493 DO l = 0, 3 494 ! 495 !-- Default surfaces 496 DO m = 1, surf_def_v(l)%ns 497 i = surf_def_v(l)%i(m) 498 j = surf_def_v(l)%j(m) 499 k = surf_def_v(l)%k(m) 500 e_p(k,j,i) = ( surf_def_v(l)%us(m) / c_0 )**2 501 ENDDO 502 ! 503 !-- Natural surfaces 504 DO m = 1, surf_lsm_v(l)%ns 505 i = surf_lsm_v(l)%i(m) 506 j = surf_lsm_v(l)%j(m) 507 k = surf_lsm_v(l)%k(m) 508 e_p(k,j,i) = ( surf_lsm_v(l)%us(m) / c_0 )**2 509 ENDDO 510 ! 511 !-- Urban surfaces 512 DO m = 1, surf_usm_v(l)%ns 513 i = surf_usm_v(l)%i(m) 514 j = surf_usm_v(l)%j(m) 515 k = surf_usm_v(l)%k(m) 516 e_p(k,j,i) = ( surf_usm_v(l)%us(m) / c_0 )**2 517 ENDDO 518 ENDDO 519 ENDIF 520 ! 521 !-- Set Neumann boundary condition for TKE at model top. Do this also 522 !-- in case of a nested run. 523 !$ACC KERNELS PRESENT(e_p) 524 e_p(nzt+1,:,:) = e_p(nzt,:,:) 525 !$ACC END KERNELS 526 ! 527 !-- Nesting case: if parent operates in RANS mode and child in LES mode, 528 !-- no TKE is transfered. This case, set Neumann conditions at lateral and 529 !-- top child boundaries. 530 !-- If not ( both either in RANS or in LES mode ), TKE boundary condition 531 !-- is treated in the nesting. 532 If ( child_domain ) THEN 533 IF ( rans_mode_parent .AND. .NOT. rans_mode ) THEN 534 535 e_p(nzt+1,:,:) = e_p(nzt,:,:) 536 IF ( bc_dirichlet_l ) e_p(:,:,nxl-1) = e_p(:,:,nxl) 537 IF ( bc_dirichlet_r ) e_p(:,:,nxr+1) = e_p(:,:,nxr) 538 IF ( bc_dirichlet_s ) e_p(:,nys-1,:) = e_p(:,nys,:) 539 IF ( bc_dirichlet_n ) e_p(:,nyn+1,:) = e_p(:,nyn,:) 540 541 ENDIF 542 ENDIF 543 ! 544 !-- At in- and outflow boundaries also set Neumann boundary conditions 545 !-- for the SGS-TKE. An exception is made for the child domain if 546 !-- both parent and child operate in RANS mode. This case no 547 !-- lateral Neumann boundary conditions will be set but Dirichlet 548 !-- conditions will be set in the nesting. 549 IF ( .NOT. child_domain .AND. .NOT. rans_mode_parent .AND. & 550 .NOT. rans_mode ) THEN 551 IF ( bc_dirichlet_s .OR. bc_radiation_s ) THEN 552 e_p(:,nys-1,:) = e_p(:,nys,:) 553 IF ( rans_tke_e ) diss_p(:,nys-1,:) = diss_p(:,nys,:) 554 ENDIF 555 IF ( bc_dirichlet_n .OR. bc_radiation_n ) THEN 556 e_p(:,nyn+1,:) = e_p(:,nyn,:) 557 IF ( rans_tke_e ) diss_p(:,nyn+1,:) = diss_p(:,nyn,:) 558 ENDIF 559 IF ( bc_dirichlet_l .OR. bc_radiation_l ) THEN 560 e_p(:,:,nxl-1) = e_p(:,:,nxl) 561 IF ( rans_tke_e ) diss_p(:,nyn+1,:) = diss_p(:,nyn,:) 562 ENDIF 563 IF ( bc_dirichlet_r .OR. bc_radiation_r ) THEN 564 e_p(:,:,nxr+1) = e_p(:,:,nxr) 565 IF ( rans_tke_e ) diss_p(:,nyn+1,:) = diss_p(:,nyn,:) 566 ENDIF 567 ENDIF 568 ENDIF 569 570 ! 571 !-- Boundary conditions for TKE dissipation rate in RANS mode. 572 IF ( rans_tke_e ) THEN 573 ! 574 !-- Use wall function within constant-flux layer 575 !-- Upward-facing surfaces 576 !-- Default surfaces 577 DO m = 1, surf_def_h(0)%ns 578 i = surf_def_h(0)%i(m) 579 j = surf_def_h(0)%j(m) 580 k = surf_def_h(0)%k(m) 581 diss_p(k,j,i) = surf_def_h(0)%us(m)**3 & 582 / ( kappa * surf_def_h(0)%z_mo(m) ) 583 ENDDO 584 ! 585 !-- Natural surfaces 586 DO m = 1, surf_lsm_h%ns 587 i = surf_lsm_h%i(m) 588 j = surf_lsm_h%j(m) 589 k = surf_lsm_h%k(m) 590 diss_p(k,j,i) = surf_lsm_h%us(m)**3 & 591 / ( kappa * surf_lsm_h%z_mo(m) ) 592 ENDDO 593 ! 594 !-- Urban surfaces 595 DO m = 1, surf_usm_h%ns 596 i = surf_usm_h%i(m) 597 j = surf_usm_h%j(m) 598 k = surf_usm_h%k(m) 599 diss_p(k,j,i) = surf_usm_h%us(m)**3 & 600 / ( kappa * surf_usm_h%z_mo(m) ) 601 ENDDO 602 ! 603 !-- Vertical surfaces 604 DO l = 0, 3 605 ! 606 !-- Default surfaces 607 DO m = 1, surf_def_v(l)%ns 608 i = surf_def_v(l)%i(m) 609 j = surf_def_v(l)%j(m) 610 k = surf_def_v(l)%k(m) 611 diss_p(k,j,i) = surf_def_v(l)%us(m)**3 & 612 / ( kappa * surf_def_v(l)%z_mo(m) ) 613 ENDDO 614 ! 615 !-- Natural surfaces 616 DO m = 1, surf_lsm_v(l)%ns 617 i = surf_lsm_v(l)%i(m) 618 j = surf_lsm_v(l)%j(m) 619 k = surf_lsm_v(l)%k(m) 620 diss_p(k,j,i) = surf_lsm_v(l)%us(m)**3 & 621 / ( kappa * surf_lsm_v(l)%z_mo(m) ) 622 ENDDO 623 ! 624 !-- Urban surfaces 625 DO m = 1, surf_usm_v(l)%ns 626 i = surf_usm_v(l)%i(m) 627 j = surf_usm_v(l)%j(m) 628 k = surf_usm_v(l)%k(m) 629 diss_p(k,j,i) = surf_usm_v(l)%us(m)**3 & 630 / ( kappa * surf_usm_v(l)%z_mo(m) ) 631 ENDDO 632 ENDDO 633 ! 634 !-- Limit change of diss to be between -90% and +100%. Also, set an absolute 635 !-- minimum value 636 DO i = nxl, nxr 637 DO j = nys, nyn 638 DO k = nzb, nzt+1 639 diss_p(k,j,i) = MAX( MIN( diss_p(k,j,i), & 640 2.0_wp * diss(k,j,i) ), & 641 0.1_wp * diss(k,j,i), & 642 0.0001_wp ) 643 ENDDO 644 ENDDO 645 ENDDO 646 647 648 diss_p(nzt+1,:,:) = diss_p(nzt,:,:) 649 650 ENDIF 651 END SUBROUTINE tcm_boundary_conds 652 373 653 !------------------------------------------------------------------------------! 374 654 ! Description: … … 926 1206 USE model_1d_mod, & 927 1207 ONLY: e1d, kh1d, km1d 928 929 USE surface_mod, &930 ONLY: get_topography_top_index_ji931 1208 932 1209 IMPLICIT NONE … … 1713 1990 ONLY: phi_m 1714 1991 1715 USE surface_mod, &1716 ONLY : surf_def_h, surf_lsm_h, surf_usm_h1717 1718 1992 IMPLICIT NONE 1719 1993 … … 2418 2692 USE bulk_cloud_model_mod, & 2419 2693 ONLY: bulk_cloud_model 2420 2421 USE surface_mod, &2422 ONLY : surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, &2423 surf_usm_v2424 2694 2425 2695 IMPLICIT NONE … … 3122 3392 ONLY: bulk_cloud_model 3123 3393 3124 USE surface_mod, &3125 ONLY : surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, &3126 surf_usm_v3127 3128 3394 IMPLICIT NONE 3129 3395 … … 3767 4033 ONLY: use_sgs_for_particles, wang_kernel 3768 4034 3769 USE surface_mod, &3770 ONLY : bc_h3771 3772 4035 IMPLICIT NONE 3773 4036 … … 3929 4192 3930 4193 ! 3931 !-- Neumann boundary condition for dissipation diss(nzb,:,:) = diss(nzb+1,:,:) 4194 !-- Neumann boundary condition for dissipation diss(nzb,:,:) = diss(nzb+1,:,:). 4195 !-- Note, bc cannot be set in tcm_boundary conditions as the dissipation 4196 !-- in LES mode is only a diagnostic quantity. 3932 4197 IF ( .NOT. rans_tke_e .AND. ( use_sgs_for_particles .OR. & 3933 4198 wang_kernel .OR. collision_turbulence ) ) THEN … … 3973 4238 USE particle_attributes, & 3974 4239 ONLY: use_sgs_for_particles, wang_kernel 3975 3976 USE surface_mod, &3977 ONLY : bc_h3978 4240 3979 4241 IMPLICIT NONE … … 4062 4324 !-- For each surface type determine start and end index (in case of elevated 4063 4325 !-- topography several up/downward facing surfaces may exist. 4326 !-- Note, bc cannot be set in tcm_boundary conditions as the dissipation 4327 !-- in LES mode is only a diagnostic quantity. 4064 4328 surf_s = bc_h(0)%start_index(j,i) 4065 4329 surf_e = bc_h(0)%end_index(j,i) … … 4322 4586 ONLY: phi_m 4323 4587 4324 USE surface_mod, &4325 ONLY : bc_h, surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, &4326 surf_usm_h, surf_usm_v4327 4328 4588 INTEGER(iwp) :: i !< loop index 4329 4589 INTEGER(iwp) :: j !< loop index -
palm/trunk/SOURCE/vertical_nesting_mod.f90
r4101 r4102 26 26 ! ----------------- 27 27 ! $Id$ 28 ! - Slightly revise setting of boundary conditions at horizontal walls, use 29 ! data-structure offset index instead of pre-calculate it for each facing 30 ! 31 ! 4101 2019-07-17 15:14:26Z gronemeier 28 32 ! remove old_dt 29 33 ! … … 2711 2715 INTEGER(iwp) :: i 2712 2716 INTEGER(iwp) :: j 2713 INTEGER(iwp) :: k 2714 INTEGER(iwp) :: kb !< variable to set respective boundary value, depends on facing. 2717 INTEGER(iwp) :: k 2715 2718 INTEGER(iwp) :: l !< running index boundary type, for up- and downward-facing walls 2716 2719 INTEGER(iwp) :: m !< running index surface elements … … 2860 2863 IF ( ibc_pt_b == 1 ) THEN 2861 2864 DO l = 0, 1 2862 ! 2863 !-- Set kb, for upward-facing surfaces value at topography top (k-1) is set, 2864 !-- for downward-facing surfaces at topography bottom (k+1). 2865 kb = MERGE( -1, 1, l == 0 ) 2866 DO m = 1, bc_h(l)%ns 2867 i = bc_h(l)%i(m) 2868 j = bc_h(l)%j(m) 2869 k = bc_h(l)%k(m) 2870 pt(k+kb,j,i) = pt(k,j,i) 2871 ENDDO 2872 ENDDO 2865 DO m = 1, bc_h(l)%ns 2866 i = bc_h(l)%i(m) 2867 j = bc_h(l)%j(m) 2868 k = bc_h(l)%k(m) 2869 pt(k+bc_h(l)%koff,j,i) = pt(k,j,i) 2870 ENDDO 2871 ENDDO 2873 2872 ENDIF 2874 2873
Note: See TracChangeset
for help on using the changeset viewer.