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.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • 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
Note: See TracChangeset for help on using the changeset viewer.