Changeset 4102 for palm/trunk/SOURCE


Ignore:
Timestamp:
Jul 17, 2019 4:00:03 PM (5 years ago)
Author:
suehring
Message:

Bugfix, set Neumann boundary conditions for the subgrid TKE at vertical walls instead of implicit Dirichlet conditions that always act as a sink term for the subgrid TKE. Therefore, add new data structure for vertical surfaces and revise the setting of the boundary grid point index space. Moreover, slightly revise setting of boundary conditions at upward- and downward facing surfaces. Finally, setting of boundary conditions for subgrid TKE and dissipation (in RANS mode) is now modularized. Update test case results.

Location:
palm/trunk/SOURCE
Files:
7 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/Makefile

    r4070 r4102  
    2525# -----------------
    2626# $Id$
     27# Remove dependency on pmc_interface for boundary_conds
     28#
     29# 4070 2019-07-03 13:51:40Z gronemeier
    2730# Add new data output modules
    2831#
     
    802805        radiation_model_mod.o
    803806boundary_conds.o: \
    804         basic_constants_and_equations_mod.o \
    805807        bulk_cloud_model_mod.o \
    806808        chemistry_model_mod.o \
    807809        mod_kinds.o \
    808810        modules.o \
    809         pmc_interface_mod.o \
    810811        salsa_mod.o \
    811812        surface_mod.o \
  • palm/trunk/SOURCE/boundary_conds.f90

    r4087 r4102  
    2525! -----------------
    2626! $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
    2734! Add comment
    2835!
     
    224231    USE arrays_3d,                                                             &
    225232        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                diss, 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,                          &
    227234               q_p, qc_p, qr_p, s, s_p, sa, sa_p, u, u_init, u_m_l, u_m_n,     &
    228235               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,  &
    229236               w, w_p, w_m_l, w_m_n, w_m_r, w_m_s
    230237
    231     USE basic_constants_and_equations_mod,                                     &
    232         ONLY:  kappa
    233 
    234238    USE bulk_cloud_model_mod,                                                  &
    235239        ONLY: bulk_cloud_model, microphysics_morrison, microphysics_seifert
     
    239243
    240244    USE control_parameters,                                                    &
    241         ONLY:  air_chemistry, bc_dirichlet_l, bc_dirichlet_n, bc_dirichlet_r,  &
     245        ONLY:  air_chemistry, bc_dirichlet_l,                                  &
    242246               bc_dirichlet_s, bc_radiation_l, bc_radiation_n, bc_radiation_r, &
    243247               bc_radiation_s, bc_pt_t_val, bc_q_t_val, bc_s_t_val,            &
    244                child_domain, constant_diffusion, coupling_mode, dt_3d,         &
     248               child_domain, coupling_mode, dt_3d,                             &
    245249               humidity, ibc_pt_b, ibc_pt_t, ibc_q_b, ibc_q_t, ibc_s_b,        &
    246250               ibc_s_t, ibc_uv_b, ibc_uv_t, intermediate_timestep_count,       &
    247251               nesting_offline, neutral, nudging, ocean_mode, passive_scalar,  &
    248                rans_mode, rans_tke_e, tsc, salsa, use_cmax
     252               tsc, salsa, use_cmax
    249253
    250254    USE grid_variables,                                                        &
     
    262266
    263267    USE pmc_interface,                                                         &
    264         ONLY : nesting_mode, rans_mode_parent
     268        ONLY : nesting_mode
    265269 
    266270    USE salsa_mod,                                                             &
     
    268272
    269273    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
    272275
    273276    USE turbulence_closure_mod,                                                &
    274         ONLY:  c_0
     277        ONLY:  tcm_boundary_conds
    275278
    276279    IMPLICIT NONE
     
    279282    INTEGER(iwp) ::  j  !< grid index y direction
    280283    INTEGER(iwp) ::  k  !< grid index z direction
    281     INTEGER(iwp) ::  kb !< variable to set respective boundary value, depends on facing.
    282284    INTEGER(iwp) ::  l  !< running index boundary type, for up- and downward-facing walls
    283285    INTEGER(iwp) ::  m  !< running index surface elements
     
    297299!-- of downward-facing surfaces.
    298300    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 )
    303301       !$OMP PARALLEL DO PRIVATE( i, j, k )
    304302       !$ACC PARALLEL LOOP PRIVATE(i, j, k) &
     
    308306          j = bc_h(l)%j(m)
    309307          k = bc_h(l)%k(m)
    310           w_p(k+kb,j,i) = 0.0_wp
     308          w_p(k+bc_h(l)%koff,j,i) = 0.0_wp
    311309       ENDDO
    312310    ENDDO
     
    341339       IF ( ibc_pt_b == 0 )  THEN
    342340          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 )
    347341             !$OMP PARALLEL DO PRIVATE( i, j, k )
    348342             DO  m = 1, bc_h(l)%ns
     
    350344                j = bc_h(l)%j(m)
    351345                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)
    353347             ENDDO
    354348          ENDDO
     
    357351       ELSEIF ( ibc_pt_b == 1 )  THEN
    358352          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 )
    363353             !$OMP PARALLEL DO PRIVATE( i, j, k )
    364354             !$ACC PARALLEL LOOP PRIVATE(i, j, k) &
     
    368358                j = bc_h(l)%j(m)
    369359                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)
    371361             ENDDO
    372362          ENDDO
     
    391381       ENDIF
    392382    ENDIF
    393 
    394 !
    395 !-- Boundary conditions for TKE.
    396 !-- Generally Neumann conditions with de/dz=0 are assumed.
    397     IF ( .NOT. constant_diffusion )  THEN
    398 
    399        IF ( .NOT. rans_mode )  THEN
    400           DO  l = 0, 1
    401 !
    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)%ns
    409                 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              ENDDO
    414           ENDDO
    415        ELSE
    416 !
    417 !--       Use wall function within constant-flux layer
    418 !--       Note, grid points listed in bc_h are not included in any calculations in RANS mode and
    419 !--       are therefore not set here.
    420 !
    421 !--       Upward-facing surfaces
    422 !--       Default surfaces
    423           DO  m = 1, surf_def_h(0)%ns
    424              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 )**2
    428           ENDDO
    429 !
    430 !--       Natural surfaces
    431           DO  m = 1, surf_lsm_h%ns
    432              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 )**2
    436           ENDDO
    437 !
    438 !--       Urban surfaces
    439           DO  m = 1, surf_usm_h%ns
    440              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 )**2
    444           ENDDO
    445 !
    446 !--       Vertical surfaces
    447           DO  l = 0, 3
    448 !
    449 !--          Default surfaces
    450              DO  m = 1, surf_def_v(l)%ns
    451                 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 )**2
    455              ENDDO
    456 !
    457 !--          Natural surfaces
    458              DO  m = 1, surf_lsm_v(l)%ns
    459                 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 )**2
    463              ENDDO
    464 !
    465 !--          Urban surfaces
    466              DO  m = 1, surf_usm_v(l)%ns
    467                 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 )**2
    471              ENDDO
    472           ENDDO
    473        ENDIF
    474 
    475        IF ( .NOT. child_domain )  THEN
    476           !$ACC KERNELS PRESENT(e_p)
    477           e_p(nzt+1,:,:) = e_p(nzt,:,:)
    478           !$ACC END KERNELS
    479 !
    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 and
    482 !--    top child boundaries.
    483 !--    If not ( both either in RANS or in LES mode ), TKE boundary condition
    484 !--    is treated in the nesting.
    485        ELSE
    486 
    487           IF ( rans_mode_parent  .AND.  .NOT. rans_mode )  THEN
    488 
    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           ENDIF
    496        ENDIF
    497     ENDIF
    498 
    499 !
    500 !-- Boundary conditions for TKE dissipation rate.
    501     IF ( rans_tke_e )  THEN
    502 !
    503 !--    Use wall function within constant-flux layer
    504 !--    Upward-facing surfaces
    505 !--    Default surfaces
    506        DO  m = 1, surf_def_h(0)%ns
    507           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        ENDDO
    513 !
    514 !--    Natural surfaces
    515        DO  m = 1, surf_lsm_h%ns
    516           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        ENDDO
    522 !
    523 !--    Urban surfaces
    524        DO  m = 1, surf_usm_h%ns
    525           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        ENDDO
    531 !
    532 !--    Vertical surfaces
    533        DO  l = 0, 3
    534 !
    535 !--       Default surfaces
    536           DO  m = 1, surf_def_v(l)%ns
    537              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           ENDDO
    543 !
    544 !--       Natural surfaces
    545           DO  m = 1, surf_lsm_v(l)%ns
    546              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           ENDDO
    552 !
    553 !--       Urban surfaces
    554           DO  m = 1, surf_usm_v(l)%ns
    555              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           ENDDO
    561        ENDDO
    562 !
    563 !--    Limit change of diss to be between -90% and +100%. Also, set an absolute
    564 !--    minimum value
    565        DO  i = nxl, nxr
    566           DO  j = nys, nyn
    567              DO  k = nzb, nzt+1
    568                 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              ENDDO
    573           ENDDO
    574        ENDDO
    575 
    576        IF ( .NOT. child_domain )  THEN
    577           diss_p(nzt+1,:,:) = diss_p(nzt,:,:)
    578        ENDIF
    579     ENDIF
    580 
    581383!
    582384!-- Boundary conditions for salinity
     
    586388!--    given.
    587389       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 )
    592390          !$OMP PARALLEL DO PRIVATE( i, j, k )
    593391          DO  m = 1, bc_h(l)%ns
     
    595393             j = bc_h(l)%j(m)
    596394             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)
    598396          ENDDO
    599397       ENDDO
     
    620418
    621419          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 )
    626420             !$OMP PARALLEL DO PRIVATE( i, j, k )
    627421             DO  m = 1, bc_h(l)%ns
     
    629423                j = bc_h(l)%j(m)
    630424                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)
    632426             ENDDO
    633427          ENDDO
     
    636430         
    637431          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 )
    642432             !$OMP PARALLEL DO PRIVATE( i, j, k )
    643433             DO  m = 1, bc_h(l)%ns
     
    645435                j = bc_h(l)%j(m)
    646436                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)
    648438             ENDDO
    649439          ENDDO
     
    662452!--       Run loop over all non-natural and natural walls. Note, in wall-datatype
    663453!--       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
    665456          !$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
    672464          ENDDO
    673465!
     
    683475!--       Run loop over all non-natural and natural walls. Note, in wall-datatype
    684476!--       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
    686479          !$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
    693487          ENDDO
    694488!
     
    711505         
    712506          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 )
    717507             !$OMP PARALLEL DO PRIVATE( i, j, k )
    718508             DO  m = 1, bc_h(l)%ns
     
    720510                j = bc_h(l)%j(m)
    721511                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)
    723513             ENDDO
    724514          ENDDO
     
    727517         
    728518          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 )
    733519             !$OMP PARALLEL DO PRIVATE( i, j, k )
    734520             DO  m = 1, bc_h(l)%ns
     
    736522                j = bc_h(l)%j(m)
    737523                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)
    739525             ENDDO
    740526          ENDDO
     
    750536       ENDIF
    751537
    752     ENDIF   
     538    ENDIF 
     539!
     540!-- Set boundary conditions for subgrid TKE and dissipation (RANS mode)
     541    CALL tcm_boundary_conds
    753542!
    754543!-- Top/bottom boundary conditions for chemical species
     
    760549!-- version) these levels are handled as a prognostic level, boundary values
    761550!-- have to be restored here.
    762 !-- For the SGS-TKE, Neumann boundary conditions are used at the inflow.
    763551    IF ( bc_dirichlet_s )  THEN
    764552       v_p(:,nys,:) = v_p(:,nys-1,:)
    765        IF ( .NOT. constant_diffusion ) e_p(:,nys-1,:) = e_p(:,nys,:)
    766     ELSEIF ( bc_dirichlet_n )  THEN
    767        IF ( .NOT. constant_diffusion ) e_p(:,nyn+1,:) = e_p(:,nyn,:)
    768553    ELSEIF ( bc_dirichlet_l ) THEN
    769554       u_p(:,:,nxl) = u_p(:,:,nxl-1)
    770        IF ( .NOT. constant_diffusion ) e_p(:,:,nxl-1) = e_p(:,:,nxl)
    771     ELSEIF ( bc_dirichlet_r )  THEN
    772        IF ( .NOT. constant_diffusion ) e_p(:,:,nxr+1) = e_p(:,:,nxr)
    773555    ENDIF
    774556
     
    777559!-- in case of nest boundaries. This must not be done in case of vertical nesting
    778560!-- 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.
    779563    IF ( nesting_mode /= 'vertical'  .OR.  nesting_offline )  THEN
    780564       IF ( bc_dirichlet_s )  THEN
     
    787571
    788572!
    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.
    790576    IF ( bc_radiation_s )  THEN
    791577       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,:)
    794578       IF ( humidity )  THEN
    795579          q_p(:,nys-1,:) = q_p(:,nys,:)
     
    806590    ELSEIF ( bc_radiation_n )  THEN
    807591       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,:)
    810592       IF ( humidity )  THEN
    811593          q_p(:,nyn+1,:) = q_p(:,nyn,:)
     
    822604    ELSEIF ( bc_radiation_l )  THEN
    823605       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)
    826606       IF ( humidity )  THEN
    827607          q_p(:,:,nxl-1) = q_p(:,:,nxl)
     
    838618    ELSEIF ( bc_radiation_r )  THEN
    839619       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)
    842620       IF ( humidity )  THEN
    843621          q_p(:,:,nxr+1) = q_p(:,:,nxr)
  • palm/trunk/SOURCE/chemistry_model_mod.f90

    r4080 r4102  
    2727! -----------------
    2828! $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
    2933! Restore accidantly removed limitation to positive values
    3034!
     
    804808    INTEGER(iwp) ::  j                            !< grid index y direction.
    805809    INTEGER(iwp) ::  k                            !< grid index z direction.
    806     INTEGER(iwp) ::  kb                           !< variable to set respective boundary value, depends on facing.
    807810    INTEGER(iwp) ::  l                            !< running index boundary type, for up- and downward-facing walls.
    808811    INTEGER(iwp) ::  m                            !< running index surface elements.
     
    839842       CASE ( 'set_bc_bottomtop' )                   
    840843!
    841 !--          Bottom boundary condtions for chemical species     
     844!--       Boundary condtions for chemical species at horizontal walls     
    842845          DO  lsp = 1, nspec                                                     
    843              IF ( ibc_cs_b == 0 )  THEN                   
     846             IF ( ibc_cs_b == 0 )  THEN
    844847                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. the
    849 !--                value at the topography bottom (k+1) is set.
    850 
    851                    kb = MERGE( -1, 1, l == 0 )
    852848                   !$OMP PARALLEL DO PRIVATE( i, j, k )
    853849                   DO  m = 1, bc_h(l)%ns
     
    855851                       j = bc_h(l)%j(m)
    856852                       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)
    858855                   ENDDO                                       
    859856                ENDDO                                       
     
    863860!--             in boundary_conds there is som extra loop over m here for passive tracer
    864861                DO  l = 0, 1
    865                    kb = MERGE( -1, 1, l == 0 )
    866862                   !$OMP PARALLEL DO PRIVATE( i, j, k )                                           
    867863                   DO m = 1, bc_h(l)%ns
     
    869865                      j = bc_h(l)%j(m)
    870866                      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)
    872869
    873870                   ENDDO
  • palm/trunk/SOURCE/salsa_mod.f90

    r4079 r4102  
    2626! -----------------
    2727! $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
    2832! Application of monotonic flux limiter for the vertical scalar advection
    2933! up to the topography top (only for the cache-optimized version at the
     
    80178021    INTEGER(iwp) ::  j    !< grid index y direction
    80188022    INTEGER(iwp) ::  k    !< grid index y direction
    8019     INTEGER(iwp) ::  kb   !< variable to set respective boundary value, depends on facing.
    80208023    INTEGER(iwp) ::  l    !< running index boundary type, for up- and downward-facing walls
    80218024    INTEGER(iwp) ::  m    !< running index surface elements
     
    80308033!--       belongs to the atmospheric grid point, therefore, set s_p at k-1
    80318034          DO  l = 0, 1
    8032 !
    8033 !--          Set kb, for upward-facing surfaces value at topography top (k-1) is
    8034 !--          set, for downward-facing surfaces at topography bottom (k+1)
    8035              kb = MERGE ( -1, 1, l == 0 )
    80368035             !$OMP PARALLEL PRIVATE( ib, ic, icc, ig, i, j, k )
    80378036             !$OMP DO
     
    80438042
    80448043                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)
    80468046                   DO  ic = 1, ncomponents_mass
    80478047                      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)
    80498050                   ENDDO
    80508051                ENDDO
    80518052                IF ( .NOT. salsa_gases_from_chem )  THEN
    80528053                   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)
    80548056                   ENDDO
    80558057                ENDIF
     
    80638065
    80648066          DO l = 0, 1
    8065 !
    8066 !--          Set kb, for upward-facing surfaces value at topography top (k-1) is
    8067 !--          set, for downward-facing surfaces at topography bottom (k+1)
    8068              kb = MERGE( -1, 1, l == 0 )
    80698067             !$OMP PARALLEL PRIVATE( ib, ic, icc, ig, i, j, k )
    80708068             !$OMP DO
     
    80768074
    80778075                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)
    80798078                   DO  ic = 1, ncomponents_mass
    80808079                      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)
    80828082                   ENDDO
    80838083                ENDDO
    80848084                IF ( .NOT. salsa_gases_from_chem ) THEN
    80858085                   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)
    80878088                   ENDDO
    80888089                ENDIF
  • palm/trunk/SOURCE/surface_mod.f90

    r3943 r4102  
    2626! -----------------
    2727! $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
    2832! Removed qsws_eb as it is no longer needed.
    2933!
     
    306310!-- are applied
    307311    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   
    314320
    315321       INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: start_index !< start index within surface data type for given (j,i)
     
    600606
    601607    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
    602609
    603610    TYPE (surf_type), DIMENSION(0:2), TARGET ::  surf_def_h  !< horizontal default surfaces (Up, Down, and Top)
     
    668675!
    669676!-- 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,      &
    672679           vertical_surfaces_exist, surf_bulk_cloud_model, surf_microphysics_morrison,             &
    673680           surf_microphysics_seifert
     
    687694! Description:
    688695! ------------
    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.
    690698!------------------------------------------------------------------------------!
    691699    SUBROUTINE init_bc
     
    696704       INTEGER(iwp) ::  j         !< loop index along y-direction
    697705       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
    721809             ENDDO
    722810          ENDDO
    723811       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
    767831                   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
    776870                   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
    778876             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) - 1
    781              start_index_h(0)         = bc_h(0)%end_index(j,i) + 1
    782 
    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) - 1
    785              start_index_h(1)         = bc_h(1)%end_index(j,i) + 1
    786877          ENDDO
    787878       ENDDO
  • palm/trunk/SOURCE/turbulence_closure_mod.f90

    r4048 r4102  
    2525! -----------------
    2626! $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
    2736! write out preprocessor directives; remove tailing whitespaces
    2837!
     
    192201!> @todo Check for random disturbances
    193202!> @note <Enter notes on the module>
    194 !------------------------------------------------------------------------------!
     203!-----------------------------------------------------------------------------!
    195204 MODULE turbulence_closure_mod
    196205
    197206
    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,          &
    201210               te_m, tend, u, v, vpt, w
    202211
    203     USE basic_constants_and_equations_mod,                                     &
     212    USE basic_constants_and_equations_mod,                                    &
    204213        ONLY:  g, kappa, lv_d_cp, lv_d_rd, rd_d_rv
    205214
    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,         &
    214232               ws_scheme_sca, current_timestep_number
    215233
    216     USE advec_ws,                                                              &
     234    USE advec_ws,                                                             &
    217235        ONLY:  advec_s_ws
    218236
    219     USE advec_s_bc_mod,                                                        &
     237    USE advec_s_bc_mod,                                                       &
    220238        ONLY:  advec_s_bc
    221239
    222     USE advec_s_pw_mod,                                                        &
     240    USE advec_s_pw_mod,                                                       &
    223241        ONLY:  advec_s_pw
    224242
    225     USE advec_s_up_mod,                                                        &
     243    USE advec_s_up_mod,                                                       &
    226244        ONLY:  advec_s_up
    227245
    228     USE cpulog,                                                                &
     246    USE cpulog,                                                               &
    229247        ONLY:  cpu_log, log_point_s
    230248
    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,    &
    233251               wall_flags_0
    234252
    235253    USE kinds
    236254
    237     USE ocean_mod,                                                             &
     255    USE ocean_mod,                                                            &
    238256        ONLY:  prho_reference
    239257
    240258    USE pegrid
    241259
    242     USE plant_canopy_model_mod,                                                &
     260    USE plant_canopy_model_mod,                                               &
    243261        ONLY:  pcm_tendency
    244262
    245     USE statistics,                                                            &
     263    USE statistics,                                                           &
    246264        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
    248276
    249277    IMPLICIT NONE
     
    271299    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  l_wall !< near-wall mixing length
    272300
    273 
     301!
     302!-- Public variables
    274303    PUBLIC c_0, rans_const_c, rans_const_sigma
    275304
     
    277306
    278307    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,                                                     &
    293323       tcm_diffusivities
    294324
    295325!
    296326!-- 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!
    297332!-- Input parameter checks to be done in check_parameters
    298333    INTERFACE tcm_check_parameters
     
    371406 CONTAINS
    372407
     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 
    373653!------------------------------------------------------------------------------!
    374654! Description:
     
    9261206    USE model_1d_mod,                                                          &
    9271207        ONLY:  e1d, kh1d, km1d
    928 
    929     USE surface_mod,                                                           &
    930         ONLY:  get_topography_top_index_ji
    9311208
    9321209    IMPLICIT NONE
     
    17131990        ONLY:  phi_m
    17141991
    1715     USE surface_mod,                                                           &
    1716         ONLY :  surf_def_h, surf_lsm_h, surf_usm_h
    1717 
    17181992    IMPLICIT NONE
    17191993
     
    24182692    USE bulk_cloud_model_mod,                                                  &
    24192693        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_v
    24242694
    24252695    IMPLICIT NONE
     
    31223392        ONLY:  bulk_cloud_model
    31233393
    3124     USE surface_mod,                                                           &
    3125         ONLY :  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h,    &
    3126                 surf_usm_v
    3127 
    31283394    IMPLICIT NONE
    31293395
     
    37674033        ONLY:  use_sgs_for_particles, wang_kernel
    37684034
    3769     USE surface_mod,                                                           &
    3770        ONLY :  bc_h
    3771 
    37724035    IMPLICIT NONE
    37734036
     
    39294192
    39304193!
    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.
    39324197    IF ( .NOT. rans_tke_e .AND. ( use_sgs_for_particles  .OR.                  &
    39334198         wang_kernel  .OR.  collision_turbulence  ) )  THEN
     
    39734238    USE particle_attributes,                                                   &
    39744239        ONLY:  use_sgs_for_particles, wang_kernel
    3975 
    3976     USE surface_mod,                                                           &
    3977        ONLY :  bc_h
    39784240
    39794241    IMPLICIT NONE
     
    40624324!--    For each surface type determine start and end index (in case of elevated
    40634325!--    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.
    40644328       surf_s = bc_h(0)%start_index(j,i)
    40654329       surf_e = bc_h(0)%end_index(j,i)
     
    43224586        ONLY:  phi_m
    43234587
    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_v
    4327 
    43284588    INTEGER(iwp) ::  i          !< loop index
    43294589    INTEGER(iwp) ::  j          !< loop index
  • palm/trunk/SOURCE/vertical_nesting_mod.f90

    r4101 r4102  
    2626! -----------------
    2727! $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
    2832! remove old_dt
    2933!
     
    27112715        INTEGER(iwp)                          ::  i
    27122716        INTEGER(iwp)                          ::  j
    2713         INTEGER(iwp)                          ::  k
    2714         INTEGER(iwp)                          ::  kb !< variable to set respective boundary value, depends on facing.
     2717        INTEGER(iwp)                          ::  k
    27152718        INTEGER(iwp)                          ::  l  !< running index boundary type, for up- and downward-facing walls
    27162719        INTEGER(iwp)                          ::  m  !< running index surface elements
     
    28602863            IF ( ibc_pt_b == 1 )  THEN
    28612864               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
    28732872            ENDIF
    28742873   
Note: See TracChangeset for help on using the changeset viewer.