Ignore:
Timestamp:
Jun 19, 2018 2:03:12 PM (3 years ago)
Author:
gronemeier
Message:

merge with branch rans

Location:
palm/trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk

  • palm/trunk/SOURCE

  • palm/trunk/SOURCE/turbulence_closure_mod.f90

    r3045 r3083  
    2525! -----------------
    2626! $Id$
     27! - set limits of diss at the end of prognostic equations
     28! - call production_e to calculate production term of diss
     29! - limit change of diss to -90% to +100%
     30! - remove factor 0.5 from diffusion_diss_ij
     31! - rename c_m into c_0, and c_h into c_4
     32! - add rans_const_c and rans_const_sigma as namelist parameters
     33! - add calculation of mixing length for profile output in case of rans_tke_e
     34! - changed format of annotations to comply with doxygen standards
     35! - calculate and save dissipation rate during rans_tke_l mode
     36! - set bc at vertical walls for e, diss, km, kh
     37! - bugfix: set l_wall = 0.0 within buildings
     38! - set l_wall at bottom and top boundary (rans-mode)
     39! - bugfix in production term for dissipation rate
     40! - bugfix in diffusion of dissipation rate
     41! - disable check for 1D model if rans_tke_e is used
     42! - bugfixes for initialization (rans-mode):
     43!    - correction of dissipation-rate formula
     44!    - calculate km based on l_wall
     45!    - initialize diss if 1D model is not used
     46!
     47! 3045 2018-05-28 07:55:41Z Giersch
    2748! Error message revised
    2849!
     
    3657! Further todo's
    3758!
    38 ! 2936 2018-03-27 14:49:27Z suehring
     59! 2936 2018-03-27 14:49:27Z gronemeier
    3960! - defined l_grid only within this module
    4061! - Moved l_wall definition from modules.f90
     
    81102!>       add OpenMP directives whereever possible
    82103!>       remove debug output variables (dummy1, dummy2, dummy3)
    83 !> @todo Move initialization of wall-mixing length from init_grid
    84104!> @todo Check for random disturbances
    85105!> @note <Enter notes on the module>
     
    146166
    147167
    148     REAL(wp) ::  c_1  = 1.44_wp    !< model constant for RANS mode
    149     REAL(wp) ::  c_2  = 1.92_wp    !< model constant for RANS mode
    150     REAL(wp) ::  c_3  = 1.44_wp    !< model constant for RANS mode
    151     REAL(wp) ::  c_h  = 0.0015_wp  !< model constant for RANS mode
    152     REAL(wp) ::  c_m               !< constant used for diffusion coefficient and dissipation (dependent on mode RANS/LES)
    153     REAL(wp) ::  c_mu = 0.09_wp    !< model constant for RANS mode
    154     REAL(wp) ::  l_max             !< maximum length scale for Blackadar mixing length
    155     REAL(wp) ::  sig_e = 1.0_wp    !< factor to calculate Ke from Km
    156     REAL(wp) ::  sig_diss = 1.3_wp !< factor to calculate K_diss from Km
    157     INTEGER(iwp) ::  surf_e        !< end index of surface elements at given i-j position
    158     INTEGER(iwp) ::  surf_s        !< start index of surface elements at given i-j position
     168    REAL(wp) ::  c_0                !< constant used for diffusion coefficient and dissipation (dependent on mode RANS/LES)
     169    REAL(wp) ::  c_1                !< model constant for RANS mode
     170    REAL(wp) ::  c_2                !< model constant for RANS mode
     171    REAL(wp) ::  c_3                !< model constant for RANS mode
     172    REAL(wp) ::  c_4                !< model constant for RANS mode
     173    REAL(wp) ::  l_max              !< maximum length scale for Blackadar mixing length
     174    REAL(wp) ::  dsig_e = 1.0_wp    !< factor to calculate Ke from Km (1/sigma_e)
     175    REAL(wp) ::  dsig_diss = 1.0_wp !< factor to calculate K_diss from Km (1/sigma_diss)
     176    INTEGER(iwp) ::  surf_e         !< end index of surface elements at given i-j position
     177    INTEGER(iwp) ::  surf_s         !< start index of surface elements at given i-j position
     178
     179    REAL(wp), DIMENSION(0:4) :: rans_const_c = &       !< model constants for RANS mode (namelist param)
     180       (/ 0.55_wp, 1.44_wp, 1.92_wp, 0.0_wp, 0.0_wp /) !> default values fit for standard-tke-e closure
     181
     182    REAL(wp), DIMENSION(2) :: rans_const_sigma = &     !< model constants for RANS mode, sigma values (sigma_e, sigma_diss) (namelist param)
     183       (/ 1.0_wp, 0.77_wp /)
    159184
    160185    REAL(wp), DIMENSION(:), ALLOCATABLE ::  l_black    !< mixing length according to Blackadar
     
    163188    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  l_wall !< near-wall mixing length
    164189
    165     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: dummy1 !< debug output variable
    166     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: dummy2 !< debug output variable
    167     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: dummy3 !< debug output variable
    168 
    169 
    170     PUBLIC c_m, dummy1, dummy2, dummy3
     190    !> @todo remove debug variables
     191    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: diss_prod1, diss_adve1, diss_diff1, &
     192                                               diss_prod2, diss_adve2, diss_diff2, &
     193                                               diss_prod3, diss_adve3, diss_diff3, &
     194                                               dummy1, dummy2, dummy3
     195
     196
     197    PUBLIC c_0, rans_const_c, rans_const_sigma
    171198
    172199!
     
    294321! ------------
    295322!> Check parameters routine for turbulence closure module.
     323!> @todo remove rans_mode from initialization namelist and rework checks
     324!>   The way it is implemented at the moment, the user has to set two variables
     325!>   so that the RANS mode is working. It would be better if only one parameter
     326!>   has to be set.
     327!>   2018-06-18, gronemeier
    296328!------------------------------------------------------------------------------!
    297329 SUBROUTINE tcm_check_parameters
     
    307339    IF ( rans_mode )  THEN
    308340
    309        c_m = 0.4_wp  !according to Detering and Etling (1985)
     341!
     342!--    Assign values to constants for RANS mode
     343       dsig_e    = 1.0_wp / rans_const_sigma(1)
     344       dsig_diss = 1.0_wp / rans_const_sigma(2)
     345
     346       c_0 = rans_const_c(0)
     347       c_1 = rans_const_c(1)
     348       c_2 = rans_const_c(2)
     349       !c_3 = rans_const_c(3)   !> @todo clarify how to switch between different models
     350       c_4 = rans_const_c(4)
    310351
    311352       SELECT CASE ( TRIM( turbulence_closure ) )
     
    316357          CASE ( 'TKE-e' )
    317358             rans_tke_e = .TRUE.
    318 
    319              IF ( INDEX( initializing_actions, 'set_1d-model_profiles' ) == 0  &
    320                   .AND.  .NOT.  nest_domain )  THEN
    321                 message_string = 'Initializing without 1D model while ' //     &
    322                                  'using TKE-e closure&' //                     &
    323                                  'is not possible at the moment!'
    324                 CALL message( 'tcm_check_parameters', 'TG0005', 1, 2, 0, 6, 0 )
    325              ENDIF
    326359
    327360          CASE DEFAULT
    328361             message_string = 'Unknown turbulence closure: ' //                &
    329362                              TRIM( turbulence_closure )
    330              CALL message( 'tcm_check_parameters', 'TG0001', 1, 2, 0, 6, 0 )
     363             CALL message( 'tcm_check_parameters', 'PA0500', 1, 2, 0, 6, 0 )
    331364
    332365       END SELECT
     366
     367       IF ( turbulent_inflow .OR. turbulent_outflow )  THEN
     368          message_string = 'turbulent inflow/outflow is not yet '//            &
     369                           'implemented for RANS mode'
     370          CALL message( 'tcm_check_parameters', 'PA0501', 1, 2, 0, 6, 0 )
     371       ENDIF
    333372
    334373       message_string = 'RANS mode is still in development! ' //               &
    335374                        '&Not all features of PALM are yet compatible '//      &
    336375                        'with RANS mode. &Use at own risk!'
    337        CALL message( 'tcm_check_parameters', 'TG0003', 0, 1, 0, 6, 0 )
     376       CALL message( 'tcm_check_parameters', 'PA0502', 0, 1, 0, 6, 0 )
    338377
    339378    ELSE
    340379
    341        c_m = 0.1_wp !according to Lilly (1967) and Deardorff (1980)
     380       c_0 = 0.1_wp !according to Lilly (1967) and Deardorff (1980)
     381
     382       dsig_e = 1.0_wp !assure to use K_m to calculate TKE instead
     383                       !of K_e which is used in RANS mode
    342384
    343385       SELECT CASE ( TRIM( turbulence_closure ) )
     
    347389
    348390          CASE DEFAULT
     391             !> @todo rework this part so that only one call of this error exists
    349392             message_string = 'Unknown turbulence closure: ' //                &
    350393                              TRIM( turbulence_closure )
    351              CALL message( 'tcm_check_parameters', 'TG0001', 1, 2, 0, 6, 0 )
     394             CALL message( 'tcm_check_parameters', 'PA0500', 1, 2, 0, 6, 0 )
    352395
    353396       END SELECT
    354 
    355     ENDIF
    356 
    357     IF ( rans_tke_e )  THEN
    358 
    359        IF ( turbulent_inflow .OR. turbulent_outflow )  THEN
    360           message_string = 'turbulent inflow/outflow is not yet '//            &
    361                            'implemented for TKE-e closure'
    362           CALL message( 'tcm_check_parameters', 'TG0002', 1, 2, 0, 6, 0 )
    363        ENDIF
    364397
    365398    ENDIF
     
    379412    IMPLICIT NONE
    380413
    381     CHARACTER (LEN=*) ::  unit     !<
    382     CHARACTER (LEN=*) ::  var      !<
    383 
    384     INTEGER(iwp) ::  i      !<
    385     INTEGER(iwp) ::  ilen   !<
    386     INTEGER(iwp) ::  k      !<
     414    CHARACTER (LEN=*) ::  unit     !< unit of output variable
     415    CHARACTER (LEN=*) ::  var      !< name of output variable
     416
     417    INTEGER(iwp) ::  i      !< index of var in data_output
     418    INTEGER(iwp) ::  ilen   !< length of var string
     419    INTEGER(iwp) ::  k      !< flag if var contains one of '_xy', '_xz' or '_yz'
    387420
    388421    SELECT CASE ( TRIM( var ) )
    389422
    390423       CASE ( 'diss' )
    391           IF ( .NOT.  rans_tke_e )  THEN
    392              message_string = 'output of "' // TRIM( var ) // '" requi' //  &
    393                       'res TKE-e closure for RANS mode.'
    394              CALL message( 'tcm_check_data_output', 'TG0101', 1, 2, 0, 6, 0 )
    395           ENDIF
    396424          unit = 'm2/s3'
    397425
    398        CASE ( 'dummy2', 'dummy3', 'dummy1' )
    399           unit = 'mixing length'
     426       CASE ( 'diss1', 'diss2',                         &                      !> @todo remove later
     427              'diss_prod1', 'diss_adve1', 'diss_diff1', &
     428              'diss_prod2', 'diss_adve2', 'diss_diff2', &
     429              'diss_prod3', 'diss_adve3', 'diss_diff3', 'dummy3'  )
     430          unit = 'debug output'
    400431
    401432       CASE ( 'kh', 'km' )
     
    420451    IMPLICIT NONE
    421452
    422     CHARACTER (LEN=*), INTENT(OUT) ::  grid_x   !<
    423     CHARACTER (LEN=*), INTENT(OUT) ::  grid_y   !<
    424     CHARACTER (LEN=*), INTENT(OUT) ::  grid_z   !<
    425     CHARACTER (LEN=*), INTENT(IN)  ::  var      !<
    426    
    427     LOGICAL, INTENT(OUT) ::  found   !<
    428    
     453    CHARACTER (LEN=*), INTENT(OUT) ::  grid_x   !< x grid of output variable
     454    CHARACTER (LEN=*), INTENT(OUT) ::  grid_y   !< y grid of output variable
     455    CHARACTER (LEN=*), INTENT(OUT) ::  grid_z   !< z grid of output variable
     456    CHARACTER (LEN=*), INTENT(IN)  ::  var      !< name of output variable
     457
     458    LOGICAL, INTENT(OUT) ::  found   !< flag if output variable is found
     459
    429460    found  = .TRUE.
    430461
     
    438469          grid_z = 'zu'
    439470
    440        CASE ( 'dummy2', 'dummy3', 'dummy1' )                                    !### remove later
     471       CASE ( 'diss1', 'diss2',                         &                       !> @todo remove later
     472              'diss_prod1', 'diss_adve1', 'diss_diff1', &
     473              'diss_prod2', 'diss_adve2', 'diss_diff2', &
     474              'diss_prod3', 'diss_adve3', 'diss_diff3', 'dummy3' )
    441475          grid_x = 'x'
    442476          grid_y = 'y'
     
    480514    IMPLICIT NONE
    481515
    482     CHARACTER (LEN=*) ::  mode       !<
    483     CHARACTER (LEN=*) ::  variable   !<
    484 
    485     INTEGER(iwp) ::  i   !<
    486     INTEGER(iwp) ::  j   !<
    487     INTEGER(iwp) ::  k   !<
     516    CHARACTER (LEN=*) ::  mode       !< flag defining mode 'allocate', 'sum' or 'average'
     517    CHARACTER (LEN=*) ::  variable   !< name of variable
     518
     519    INTEGER(iwp) ::  i   !< loop index
     520    INTEGER(iwp) ::  j   !< loop index
     521    INTEGER(iwp) ::  k   !< loop index
    488522
    489523    IF ( mode == 'allocate' )  THEN
     
    616650    IMPLICIT NONE
    617651
    618     CHARACTER (LEN=*) ::  grid       !<
    619     CHARACTER (LEN=*) ::  mode       !<
    620     CHARACTER (LEN=*) ::  variable   !<
    621 
    622     INTEGER(iwp) ::  av   !<
    623     INTEGER(iwp) ::  i    !<
    624     INTEGER(iwp) ::  j    !<
    625     INTEGER(iwp) ::  k    !<
    626     INTEGER(iwp) ::  nzb_do   !<
    627     INTEGER(iwp) ::  nzt_do   !<
    628 
    629     LOGICAL ::  found   !<
     652    CHARACTER (LEN=*) ::  grid       !< name of vertical grid
     653    CHARACTER (LEN=*) ::  mode       !< either 'xy', 'xz' or 'yz'
     654    CHARACTER (LEN=*) ::  variable   !< name of variable
     655
     656    INTEGER(iwp) ::  av   !< flag for (non-)average output
     657    INTEGER(iwp) ::  i    !< loop index
     658    INTEGER(iwp) ::  j    !< loop index
     659    INTEGER(iwp) ::  k    !< loop index
     660    INTEGER(iwp) ::  nzb_do   !< vertical output index (bottom)
     661    INTEGER(iwp) ::  nzt_do   !< vertical output index (top)
     662
     663    LOGICAL ::  found   !< flag if output variable is found
    630664    LOGICAL ::  two_d   !< flag parameter that indicates 2D variables (horizontal cross sections)
    631665
     
    737771    IMPLICIT NONE
    738772
    739     CHARACTER (LEN=*) ::  variable   !<
    740 
    741     INTEGER(iwp) ::  av     !<
    742     INTEGER(iwp) ::  i      !<
    743     INTEGER(iwp) ::  j      !<
    744     INTEGER(iwp) ::  k      !<
     773    CHARACTER (LEN=*) ::  variable   !< name of variable
     774
     775    INTEGER(iwp) ::  av     !< flag for (non-)average output
     776    INTEGER(iwp) ::  i      !< loop index
     777    INTEGER(iwp) ::  j      !< loop index
     778    INTEGER(iwp) ::  k      !< loop index
    745779    INTEGER(iwp) ::  nzb_do !< lower limit of the data output (usually 0)
    746780    INTEGER(iwp) ::  nzt_do !< vertical upper limit of the data output (usually nz_do3d)
    747781
    748     LOGICAL ::  found   !<
     782    LOGICAL ::  found   !< flag if output variable is found
    749783
    750784    REAL(wp) ::  fill_value = -999.0_wp  !< value for the _FillValue attribute
     
    829863          ENDIF
    830864
    831        CASE ( 'dummy1' )                                                        !### remove later
     865       CASE ( 'dummy3' )                                                        !> @todo remove later
     866          IF ( av == 0 )  THEN
     867             DO  i = nxl, nxr
     868                DO  j = nys, nyn
     869                   DO  k = nzb_do, nzt_do
     870                      local_pf(i,j,k) = dummy3(k,j,i)
     871                   ENDDO
     872                ENDDO
     873             ENDDO
     874          ENDIF
     875
     876       CASE ( 'diss1' )                                                         !> @todo remove later
    832877          IF ( av == 0 )  THEN
    833878             DO  i = nxl, nxr
     
    840885          ENDIF
    841886
    842        CASE ( 'dummy2' )                                                        !### remove later
     887       CASE ( 'diss2' )                                                         !> @todo remove later
    843888          IF ( av == 0 )  THEN
    844889             DO  i = nxl, nxr
     
    851896          ENDIF
    852897
    853        CASE ( 'dummy3' )                                                        !### remove later
     898       CASE ( 'diss_prod1' )                                                    !> @todo remove later
    854899          IF ( av == 0 )  THEN
    855900             DO  i = nxl, nxr
    856901                DO  j = nys, nyn
    857902                   DO  k = nzb_do, nzt_do
    858                       local_pf(i,j,k) = dummy3(k,j,i)
     903                      local_pf(i,j,k) = diss_prod1(k,j,i)
    859904                   ENDDO
    860905                ENDDO
     
    862907          ENDIF
    863908
     909       CASE ( 'diss_adve1' )                                                    !> @todo remove later
     910          IF ( av == 0 )  THEN
     911             DO  i = nxl, nxr
     912                DO  j = nys, nyn
     913                   DO  k = nzb_do, nzt_do
     914                      local_pf(i,j,k) = diss_adve1(k,j,i)
     915                   ENDDO
     916                ENDDO
     917             ENDDO
     918          ENDIF
     919
     920       CASE ( 'diss_diff1' )                                                    !> @todo remove later
     921          IF ( av == 0 )  THEN
     922             DO  i = nxl, nxr
     923                DO  j = nys, nyn
     924                   DO  k = nzb_do, nzt_do
     925                      local_pf(i,j,k) = diss_diff1(k,j,i)
     926                   ENDDO
     927                ENDDO
     928             ENDDO
     929          ENDIF
     930
     931       CASE ( 'diss_prod2' )                                                    !> @todo remove later
     932          IF ( av == 0 )  THEN
     933             DO  i = nxl, nxr
     934                DO  j = nys, nyn
     935                   DO  k = nzb_do, nzt_do
     936                      local_pf(i,j,k) = diss_prod2(k,j,i)
     937                   ENDDO
     938                ENDDO
     939             ENDDO
     940          ENDIF
     941
     942       CASE ( 'diss_adve2' )                                                    !> @todo remove later
     943          IF ( av == 0 )  THEN
     944             DO  i = nxl, nxr
     945                DO  j = nys, nyn
     946                   DO  k = nzb_do, nzt_do
     947                      local_pf(i,j,k) = diss_adve2(k,j,i)
     948                   ENDDO
     949                ENDDO
     950             ENDDO
     951          ENDIF
     952
     953       CASE ( 'diss_diff2' )                                                    !> @todo remove later
     954          IF ( av == 0 )  THEN
     955             DO  i = nxl, nxr
     956                DO  j = nys, nyn
     957                   DO  k = nzb_do, nzt_do
     958                      local_pf(i,j,k) = diss_diff2(k,j,i)
     959                   ENDDO
     960                ENDDO
     961             ENDDO
     962          ENDIF
     963
     964       CASE ( 'diss_prod3' )                                                    !> @todo remove later
     965          IF ( av == 0 )  THEN
     966             DO  i = nxl, nxr
     967                DO  j = nys, nyn
     968                   DO  k = nzb_do, nzt_do
     969                      local_pf(i,j,k) = diss_prod3(k,j,i)
     970                   ENDDO
     971                ENDDO
     972             ENDDO
     973          ENDIF
     974
     975       CASE ( 'diss_adve3' )                                                    !> @todo remove later
     976          IF ( av == 0 )  THEN
     977             DO  i = nxl, nxr
     978                DO  j = nys, nyn
     979                   DO  k = nzb_do, nzt_do
     980                      local_pf(i,j,k) = diss_adve3(k,j,i)
     981                   ENDDO
     982                ENDDO
     983             ENDDO
     984          ENDIF
     985
     986       CASE ( 'diss_diff3' )                                                    !> @todo remove later
     987          IF ( av == 0 )  THEN
     988             DO  i = nxl, nxr
     989                DO  j = nys, nyn
     990                   DO  k = nzb_do, nzt_do
     991                      local_pf(i,j,k) = diss_diff3(k,j,i)
     992                   ENDDO
     993                ENDDO
     994             ENDDO
     995          ENDIF
     996         
    864997       CASE DEFAULT
    865998          found = .FALSE.
     
    8911024    ALLOCATE( km(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    8921025
    893     ALLOCATE( dummy1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )                           !### remove later
     1026    ALLOCATE( dummy1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )                           !> @todo remove later
    8941027    ALLOCATE( dummy2(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    8951028    ALLOCATE( dummy3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    896 
    897     IF ( rans_mode )  ALLOCATE( l_black(nzb:nzt+1) )
     1029    ALLOCATE( diss_adve1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     1030    ALLOCATE( diss_adve2(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     1031    ALLOCATE( diss_adve3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     1032    ALLOCATE( diss_prod1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     1033    ALLOCATE( diss_prod2(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     1034    ALLOCATE( diss_prod3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     1035    ALLOCATE( diss_diff1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     1036    ALLOCATE( diss_diff2(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     1037    ALLOCATE( diss_diff3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     1038    dummy1 = 0.0_wp
     1039    dummy2 = 0.0_wp
     1040    dummy3 = 0.0_wp
     1041    diss_adve1 = 0.0_wp
     1042    diss_adve2 = 0.0_wp
     1043    diss_adve3 = 0.0_wp
     1044    diss_prod1 = 0.0_wp
     1045    diss_prod2 = 0.0_wp
     1046    diss_prod3 = 0.0_wp
     1047    diss_diff1 = 0.0_wp
     1048    diss_diff2 = 0.0_wp
     1049    diss_diff3 = 0.0_wp
    8981050
    8991051#if defined( __nopointer )
     
    9121064!-- they do not necessarily need to be transferred, which is attributed to
    9131065!-- the design of the model coupler which allocates memory for each variable.
    914     IF ( rans_tke_e  .OR.  use_sgs_for_particles  .OR.  wang_kernel  .OR.      &
     1066    IF ( rans_mode  .OR.  use_sgs_for_particles  .OR.  wang_kernel  .OR.       &
    9151067         collision_turbulence  .OR.  nested_run )  THEN
    9161068#if defined( __nopointer )
     
    9341086    e  => e_1;   e_p  => e_2;   te_m  => e_3
    9351087
    936     IF ( rans_tke_e  .OR.  use_sgs_for_particles  .OR.     &
     1088    IF ( rans_mode  .OR.  use_sgs_for_particles  .OR.     &
    9371089         wang_kernel  .OR.  collision_turbulence  .OR.  nested_run )  THEN
    9381090       diss => diss_1
     
    9671119    INTEGER(iwp) :: j            !< loop index
    9681120    INTEGER(iwp) :: k            !< loop index
    969     INTEGER(iwp) :: nz_s_shift   !<
    970     INTEGER(iwp) :: nz_s_shift_l !<
     1121    INTEGER(iwp) :: nz_s_shift   !< lower shift index for scalars
     1122    INTEGER(iwp) :: nz_s_shift_l !< local lower shift index in case of turbulent inflow
    9711123
    9721124!
    9731125!-- Initialize mixing length
    9741126    CALL tcm_init_mixing_length
     1127    dummy3 = l_wall                 !> @todo remove later
    9751128
    9761129!
     
    9951148
    9961149          IF ( rans_tke_e )  THEN
    997              IF ( dissipation_1d == 'prognostic' )  THEN    !### Why must this be checked?
    998                 DO  i = nxlg, nxrg                          !### Should 'diss' not always
    999                    DO  j = nysg, nyng                       !### be prognostic in case rans_tke_e?
     1150             IF ( dissipation_1d == 'prognostic' )  THEN    !> @query Why must this be checked?
     1151                DO  i = nxlg, nxrg                          ! Should 'diss' not always
     1152                   DO  j = nysg, nyng                       ! be prognostic in case rans_tke_e?
    10001153                      diss(:,j,i) = diss1d
    10011154                   ENDDO
     
    10051158                   DO  j = nysg, nyng
    10061159                      DO  k = nzb+1, nzt
    1007                          diss(k,j,i) = e(k,j,i) * SQRT( e(k,j,i) ) / l1d(k)
     1160                         diss(k,j,i) = c_0**4 * e(k,j,i)**2 / km1d(k)
    10081161                      ENDDO
    10091162                   ENDDO
     
    10161169
    10171170          IF ( constant_diffusion )  THEN
    1018              km   = km_constant
    1019              kh   = km / prandtl_number
    1020              e    = 0.0_wp
     1171             km = km_constant
     1172             kh = km / prandtl_number
     1173             e  = 0.0_wp
    10211174          ELSEIF ( e_init > 0.0_wp )  THEN
    1022              DO  k = nzb+1, nzt
    1023                 km(k,:,:) = c_m * l_grid(k) * SQRT( e_init )
     1175             DO  i = nxlg, nxrg
     1176                DO  j = nysg, nyng
     1177                   DO  k = nzb+1, nzt
     1178                      km(k,j,i) = c_0 * l_wall(k,j,i) * SQRT( e_init )
     1179                   ENDDO
     1180                ENDDO
    10241181             ENDDO
    10251182             km(nzb,:,:)   = km(nzb+1,:,:)
    10261183             km(nzt+1,:,:) = km(nzt,:,:)
    1027              kh   = km / prandtl_number
    1028              e    = e_init
     1184             kh = km / prandtl_number
     1185             e  = e_init
    10291186          ELSE
    10301187             IF ( .NOT. ocean )  THEN
     
    10401197          ENDIF
    10411198
     1199          IF ( rans_tke_e )  THEN
     1200             DO  i = nxlg, nxrg
     1201                DO  j = nysg, nyng
     1202                   DO  k = nzb+1, nzt
     1203                      diss(k,j,i) = c_0**4 * e(k,j,i)**2 / km(k,j,i)
     1204                   ENDDO
     1205                ENDDO
     1206             ENDDO
     1207             diss(nzb,:,:) = diss(nzb+1,:,:)
     1208             diss(nzt+1,:,:) = diss(nzt,:,:)
     1209          ENDIF
     1210
    10421211       ENDIF
    10431212!
     
    10731242             ENDDO
    10741243          ENDDO
     1244          IF ( rans_tke_e )  THEN
     1245             DO  i = nxlg, nxrg
     1246                DO  j = nysg, nyng
     1247                   nz_s_shift = get_topography_top_index_ji( j, i, 's' )
     1248
     1249                   diss(nz_s_shift:nzt+1,j,i) = diss(0:nzt+1-nz_s_shift,j,i)
     1250                ENDDO
     1251             ENDDO
     1252          ENDIF
    10751253       ENDIF
    10761254
     
    10841262!--       boundary and adjust mean inflow profiles
    10851263          IF ( complex_terrain )  THEN
    1086              IF ( nxlg <= 0 .AND. nxrg >= 0 .AND. nysg <= 0 .AND. nyng >= 0 )  THEN
     1264             IF ( nxlg <= 0 .AND. nxrg >= 0 .AND.  &
     1265                  nysg <= 0 .AND. nyng >= 0        )  THEN
    10871266                nz_s_shift_l = get_topography_top_index_ji( 0, 0, 's' )
    10881267             ELSE
     
    10951274             nz_s_shift = nz_s_shift_l
    10961275#endif
    1097              mean_inflow_profiles(nz_s_shift:nzt+1,5) = hom_sum(0:nzt+1-nz_s_shift,8,0)  ! e
     1276             mean_inflow_profiles(nz_s_shift:nzt+1,5) =  &
     1277                hom_sum(0:nzt+1-nz_s_shift,8,0)  ! e
    10981278          ENDIF
    10991279!
     
    11141294!
    11151295!--       Inside buildings set TKE back to zero.
    1116 !--       Other scalars (km, kh, diss, ...) are ignored at present,
     1296!--       Other scalars (km, kh,...) are ignored at present,
    11171297!--       maybe revise later.
    11181298          DO  i = nxlg, nxrg
     
    11211301                   e(k,j,i)     = MERGE( e(k,j,i), 0.0_wp,                     &
    11221302                                         BTEST( wall_flags_0(k,j,i), 0 ) )
    1123                    te_m(k,j,i)  = MERGE( te_m(k,j,i), 0.0_wp,                  &
    1124                                          BTEST( wall_flags_0(k,j,i), 0 ) )
    11251303                ENDDO
    11261304             ENDDO
    11271305          ENDDO
    11281306
     1307          IF ( rans_tke_e )  THEN
     1308             DO  i = nxlg, nxrg
     1309                DO  j = nysg, nyng
     1310                   DO  k = nzb, nzt
     1311                      diss(k,j,i)    = MERGE( diss(k,j,i), 0.0_wp,             &
     1312                                              BTEST( wall_flags_0(k,j,i), 0 ) )
     1313                   ENDDO
     1314                ENDDO
     1315             ENDDO
     1316          ENDIF
    11291317       ENDIF
    11301318!
     
    11341322!
    11351323!--    Allthough tendency arrays are set in prognostic_equations, they have
    1136 !--    have to be predefined here because they are used (but multiplied with 0)
    1137 !--    there before they are set.
     1324!--    to be predefined here because there they are used (but multiplied with 0)
     1325!--    before they are set.
    11381326       te_m = 0.0_wp
     1327
     1328       IF ( rans_tke_e )  THEN
     1329          diss_p = diss
     1330          tdiss_m = 0.0_wp
     1331       ENDIF
    11391332
    11401333    ENDIF
     
    12151408
    12161409       DO  k = 1, nzt
    1217           IF ( l_grid(k) > 1.5_wp * dx * wall_adjustment_factor .OR.           &
     1410          IF ( l_grid(k) > 1.5_wp * dx * wall_adjustment_factor .OR.            &
    12181411               l_grid(k) > 1.5_wp * dy * wall_adjustment_factor )  THEN
    1219              WRITE( message_string, * ) 'grid anisotropy exceeds ',            &
    1220                                         'threshold given by only local',       &
    1221                                         ' horizontal reduction of near_wall ', &
    1222                                         'mixing length l_wall',                &
    1223                                         ' starting from height level k = ', k, &
     1412             WRITE( message_string, * ) 'grid anisotropy exceeds ',             &
     1413                                        'threshold given by only local',        &
     1414                                        ' &horizontal reduction of near_wall ', &
     1415                                        'mixing length l_wall',                 &
     1416                                        ' &starting from height level k = ', k, &
    12241417                                        '.'
    12251418             CALL message( 'init_grid', 'PA0202', 0, 1, 0, 6, 0 )
     
    13201513!
    13211514!-- Initialize the mixing length in case of a RANS simulation
     1515       ALLOCATE( l_black(nzb:nzt+1) )
    13221516
    13231517!
    13241518!--    Calculate mixing length according to Blackadar (1962)
    13251519       IF ( f /= 0.0_wp )  THEN
    1326           l_max = 2.7E-4 * SQRT( ug(nzt+1)**2 + vg(nzt+1)**2 ) /               &
    1327                   ABS( f ) + 1E-10_wp
     1520          l_max = 2.7E-4_wp * SQRT( ug(nzt+1)**2 + vg(nzt+1)**2 ) /            &
     1521                  ABS( f ) + 1.0E-10_wp
    13281522       ELSE
    13291523          l_max = 30.0_wp
     
    13381532!
    13391533!--    Gather topography information of whole domain
    1340        !## TODO: reduce amount of data sent by MPI call
    1341        !##  By now, a whole global 3D-array is sent and received with
    1342        !##  MPI_ALLREDUCE although most of the array is 0. This can be
    1343        !##  drastically reduced if only the local subarray is sent and stored
    1344        !##  in a global array. For that, an MPI data type or subarray must be
    1345        !##  defined.
    1346        !##  2018-03-19, gronemeier
     1534       !> @todo reduce amount of data sent by MPI call
     1535       !>   By now, a whole global 3D-array is sent and received with
     1536       !>   MPI_ALLREDUCE although most of the array is 0. This can be
     1537       !>   drastically reduced if only the local subarray is sent and stored
     1538       !>   in a global array. For that, an MPI data type or subarray must be
     1539       !>   defined.
     1540       !>   2018-03-19, gronemeier
    13471541       ALLOCATE( wall_flags_0_global(nzb:nzt+1,0:ny,0:nx) )
    13481542
     
    13731567          ENDDO
    13741568       ENDDO
     1569
     1570       l_wall(nzb,:,:) = l_black(nzb)
     1571       l_wall(nzt+1,:,:) = l_black(nzt+1)
    13751572!
    13761573!--    Limit mixing length to either nearest wall or Blackadar mixing length.
     
    14151612          IF ( rad_k_b /= 0 .OR. rad_k_t /= 0 )  THEN
    14161613
    1417              !## NOTE: shape of vicinity is larger in z direction
    1418              !##  Shape of vicinity is two grid points larger than actual search
    1419              !##  radius in vertical direction. The first and last grid point is
    1420              !##   always set to 1 to asure correct detection of topography. See
    1421              !##  function "shortest_distance" for details.
    1422              !##  2018-03-16, gronemeier
     1614             !> @note shape of vicinity is larger in z direction
     1615             !>   Shape of vicinity is two grid points larger than actual search
     1616             !>   radius in vertical direction. The first and last grid point is
     1617             !>   always set to 1 to asure correct detection of topography. See
     1618             !>   function "shortest_distance" for details.
     1619             !>   2018-03-16, gronemeier
    14231620             ALLOCATE( vicinity(-rad_k-1:rad_k+1,-rad_j:rad_j,-rad_i:rad_i) )
    14241621             ALLOCATE( vic_yz(0:rad_k+1,0:rad_j) )
     
    16011798                   ELSE  !Check if (i,j,k) belongs to atmosphere
    16021799
    1603                       l_wall(k,j,i) = -999.0
     1800                      l_wall(k,j,i) = l_black(k)
    16041801
    16051802                   ENDIF
     
    16301827!> (pos_i/jj/kk), where (jj/kk) is the position of the maximum of 'array'
    16311828!> closest to the origin (0/0) of 'array'.
     1829!> @todo this part of PALM does not reproduce the same results for optimized
     1830!>   and debug options for the compiler. This should be fixed
    16321831!------------------------------------------------------------------------------!
    16331832    REAL FUNCTION shortest_distance( array, orientation, pos_i )
     
    16591858       IF ( orientation ) THEN  !if array is oriented upwards
    16601859          DO  jj = 0, rad_j
    1661              shortest_distance = MIN( shortest_distance,                        &
    1662                                       SQRT( MAX(pos_i*dx-0.5*dx,0.0)**2         &
    1663                                           + MAX(jj*dy-0.5*dy,0.0)**2            &
    1664                                           + MAX(zw(loc_k(jj)+k-1)-zu(k),0.0)**2 &
    1665                                           )                                     &
    1666                                     )
     1860             shortest_distance =                                               &
     1861                MIN( shortest_distance,                                        &
     1862                     SQRT( MAX(REAL(pos_i, KIND=wp)*dx-0.5_wp*dx, 0.0_wp)**2   &
     1863                         + MAX(REAL(jj, KIND=wp)*dy-0.5_wp*dy, 0.0_wp)**2      &
     1864                         + MAX(zw(loc_k(jj)+k-1)-zu(k), 0.0_wp)**2             &
     1865                         )                                                     &
     1866                   )
    16671867          ENDDO
    16681868       ELSE  !if array is oriented downwards
    1669           !## NOTE: MAX within zw required to circumvent error at domain border
    1670           !##  At the domain border, if non-cyclic boundary is present, the
    1671           !##  index for zw could be -1, which will be errorneous (zw(-1) does
    1672           !##  not exist). The MAX function limits the index to be at least 0.
     1869          !> @note MAX within zw required to circumvent error at domain border
     1870          !>   At the domain border, if non-cyclic boundary is present, the
     1871          !>   index for zw could be -1, which will be errorneous (zw(-1) does
     1872          !>   not exist). The MAX function limits the index to be at least 0.
    16731873          DO  jj = 0, rad_j
    1674              shortest_distance = MIN( shortest_distance,                       &
    1675                                       SQRT( MAX(pos_i*dx-0.5*dx,0.0)**2        &
    1676                                           + MAX(jj*dy-0.5*dy,0.0)**2           &
    1677                                           + MAX(zu(k)-zw(MAX(k-loc_k(jj),      &
    1678                                                              0_iwp)),          &
    1679                                                 0.0)**2                        &
    1680                                           )                                    &
    1681                                     )
     1874             shortest_distance =                                               &
     1875                MIN( shortest_distance,                                        &
     1876                     SQRT( MAX(REAL(pos_i, KIND=wp)*dx-0.5_wp*dx, 0.0_wp)**2   &
     1877                         + MAX(REAL(jj, KIND=wp)*dy-0.5_wp*dy, 0.0_wp)**2      &
     1878                         + MAX(zu(k)-zw(MAX(k-loc_k(jj),0_iwp)), 0.0_wp)**2    &
     1879                         )                                                     &
     1880                   )
    16821881          ENDDO
    16831882       ENDIF
    1684 
     1883       
    16851884    END FUNCTION
    16861885
     
    20042203!
    20052204!--    Use special boundary condition in case of TKE-e closure
     2205       !> @todo do the same for usm and lsm surfaces
     2206       !>   2018-06-05, gronemeier
    20062207       IF ( rans_tke_e )  THEN
    20072208          DO  i = nxl, nxr
     
    20112212                DO  m = surf_s, surf_e
    20122213                   k = surf_def_h(0)%k(m)
    2013                    e_p(k,j,i) = surf_def_h(0)%us(m)**2 / c_m**2
     2214                   e_p(k,j,i) = surf_def_h(0)%us(m)**2 / c_0**2
    20142215                ENDDO
    20152216             ENDDO
     
    20922293             DO  k = nzb+1, nzt
    20932294!                tend(k,j,i) = tend(k,j,i) + c_1 * diss(k,j,i) / ( e(k,j,i) + 1.0E-20_wp ) * produc(k)
    2094                 tend(k,j,i) = tend(k,j,i) + c_1 * c_mu * f / c_h               &  !### needs revision
     2295                tend(k,j,i) = tend(k,j,i) + c_1 * c_0**4 * f / c_4               &  !> @todo needs revision
    20952296                      / surf_def_h(0)%us(surf_def_h(0)%start_index(j,i))       &
    20962297                      * SQRT(e(k,j,i)) * produc(k,j,i)
     
    21032304!
    21042305!--    Additional sink term for flows through plant canopies
    2105 !        IF ( plant_canopy )  CALL pcm_tendency( ? )                            !### what to do with this?
    2106 
    2107 !        CALL user_actions( 'diss-tendency' )                                   !### not yet implemented
     2306!        IF ( plant_canopy )  CALL pcm_tendency( ? )                            !> @query what to do with this?
     2307
     2308!        CALL user_actions( 'diss-tendency' )                                   !> @todo not yet implemented
    21082309
    21092310!
     
    21812382    USE arrays_3d,                                                             &
    21822383        ONLY:  ddzu, diss_l_diss, diss_l_e, diss_s_diss, diss_s_e,             &
    2183                flux_l_diss, flux_l_e, flux_s_diss, flux_s_e
     2384               flux_l_diss, flux_l_e, flux_s_diss, flux_s_e,&
     2385               u_p,v_p,w_p
    21842386
    21852387    USE control_parameters,                                                    &
    21862388        ONLY:  f, tsc
     2389
     2390    USE grid_variables,                                                        &
     2391        ONLY:  dx, dy
    21872392
    21882393    USE surface_mod,                                                           &
     
    21902395                surf_usm_v
    21912396
     2397    use indices, only: nx, ny
     2398
    21922399    IMPLICIT NONE
    21932400
    21942401    INTEGER(iwp) ::  i       !< loop index x direction
    2195     INTEGER(iwp) ::  i_omp   !<
     2402    INTEGER(iwp) ::  i_omp   !< first loop index of i-loop in prognostic_equations
    21962403    INTEGER(iwp) ::  j       !< loop index y direction
    21972404    INTEGER(iwp) ::  k       !< loop index z direction
     2405    INTEGER(iwp) ::  l       !< loop index
    21982406    INTEGER(iwp) ::  m       !< loop index
    21992407    INTEGER(iwp) ::  surf_e  !< end index of surface elements at given i-j position
    22002408    INTEGER(iwp) ::  surf_s  !< start index of surface elements at given i-j position
    2201     INTEGER(iwp) ::  tn      !<
    2202 
    2203     REAL(wp), DIMENSION(nzb:nzt+1) :: advec  !< advection term of TKE tendency
    2204     REAL(wp), DIMENSION(nzb:nzt+1) :: produc !< production term of TKE tendency
     2409    INTEGER(iwp) ::  tn      !< task number of openmp task
     2410
     2411    INTEGER(iwp) :: pis = 32 !< debug variable, print from i=pis                !> @todo remove later
     2412    INTEGER(iwp) :: pie = 32 !< debug variable, print until i=pie               !> @todo remove later
     2413    INTEGER(iwp) :: pjs = 26 !< debug variable, print from j=pjs                !> @todo remove later
     2414    INTEGER(iwp) :: pje = 26 !< debug variable, print until j=pje               !> @todo remove later
     2415    INTEGER(iwp) :: pkb = 1  !< debug variable, print from k=pkb                !> @todo remove later
     2416    INTEGER(iwp) :: pkt = 7  !< debug variable, print until k=pkt               !> @todo remove later
     2417
     2418    REAL(wp), DIMENSION(nzb:nzt+1) :: dum_adv   !< debug variable               !> @todo remove later
     2419    REAL(wp), DIMENSION(nzb:nzt+1) :: dum_pro   !< debug variable               !> @todo remove later
     2420    REAL(wp), DIMENSION(nzb:nzt+1) :: dum_dif   !< debug variable               !> @todo remove later
     2421
     24225555 FORMAT(A,7(1X,E12.5))   !> @todo remove later
    22052423
    22062424!
     
    22242442       ENDIF
    22252443
    2226        advec(:) = tend(:,j,i)
    2227 
    2228        CALL production_e( i, j )
    2229 
    2230        produc(:) = tend(:,j,i) - advec(:)
     2444       dum_adv = tend(:,j,i)                                                    !> @todo remove later
     2445
     2446       CALL production_e( i, j, .FALSE. )
     2447
     2448       dum_pro = tend(:,j,i) - dum_adv                                          !> @todo remove later
    22312449
    22322450       IF ( .NOT. humidity )  THEN
     
    22392457          CALL diffusion_e( i, j, vpt, pt_reference )
    22402458       ENDIF
     2459
     2460       dum_dif = tend(:,j,i) - dum_adv - dum_pro                                !> @todo remove later
    22412461
    22422462!
     
    22682488          DO  m = surf_s, surf_e
    22692489             k = surf_def_h(0)%k(m)
    2270              e_p(k,j,i) = surf_def_h(0)%us(m)**2 / c_m**2
     2490             e_p(k,j,i) = ( surf_def_h(0)%us(m) / c_0 )**2
     2491          ENDDO
     2492
     2493          DO  l = 0, 3
     2494             surf_s = surf_def_v(l)%start_index(j,i)
     2495             surf_e = surf_def_v(l)%end_index(j,i)
     2496             DO  m = surf_s, surf_e
     2497                k = surf_def_v(l)%k(m)
     2498                e_p(k,j,i) = ( surf_def_v(l)%us(m) / c_0 )**2
     2499             ENDDO
    22712500          ENDDO
    22722501       ENDIF
     
    22872516          ENDIF
    22882517       ENDIF
     2518
     2519!        if ( i >= pis .and. i <= pie .and. j >= pjs .and. j <= pje ) then        !> @todo remove later
     2520!           WRITE(9, *) '------'
     2521!           WRITE(9, '(A,F8.3,1X,F8.3,1X,I2)') 't, dt, int_ts:', simulated_time, dt_3d, intermediate_timestep_count
     2522!           WRITE(9, *) 'i:', i
     2523!           WRITE(9, *) 'j:', j
     2524!           WRITE(9, *) 'k:', pkb, ' - ', pkt
     2525!           WRITE(9, *) '---'
     2526!           WRITE(9, *) 'e:'
     2527!           WRITE(9, 5555) 'adv :', dum_adv(pkb:pkt)
     2528!           WRITE(9, 5555) 'pro :', dum_pro(pkb:pkt)
     2529!           WRITE(9, 5555) 'dif :', dum_dif(pkb:pkt)
     2530!           WRITE(9, 5555) 'tend:', tend(pkb:pkt,j,i)
     2531!           WRITE(9, 5555) 'e_p :', e_p(pkb:pkt,j,i)
     2532!           WRITE(9, 5555) 'e   :', e(pkb:pkt,j,i)
     2533!           FLUSH(9)
     2534!        endif
    22892535
    22902536    ENDIF   ! TKE equation
     
    23092555       ENDIF
    23102556
     2557       IF ( intermediate_timestep_count == 1 )  diss_adve1(:,j,i) = tend(:,j,i) !> @todo remove later
     2558       IF ( intermediate_timestep_count == 2 )  diss_adve2(:,j,i) = tend(:,j,i)
     2559       IF ( intermediate_timestep_count == 3 )  diss_adve3(:,j,i) = tend(:,j,i)
     2560
    23112561!
    23122562!--    Production of TKE dissipation rate
    2313        DO  k = nzb+1, nzt
    2314 !              tend(k,j,i) = tend(k,j,i) + c_1 * diss(k,j,i) / ( e(k,j,i) + 1.0E-20_wp ) * produc(k)
    2315           tend(k,j,i) = tend(k,j,i) + c_1 * c_mu * f / c_h                     &  !### needs revision
    2316                       / surf_def_h(0)%us(surf_def_h(0)%start_index(j,i))       &
    2317                       * SQRT(e(k,j,i)) * produc(k)
    2318        ENDDO
    2319 
     2563       CALL production_e( i, j, .TRUE. )
     2564
     2565       IF ( intermediate_timestep_count == 1 )  diss_prod1(:,j,i) = tend(:,j,i) - diss_adve1(:,j,i) !> @todo remove later
     2566       IF ( intermediate_timestep_count == 2 )  diss_prod2(:,j,i) = tend(:,j,i) - diss_adve2(:,j,i)
     2567       IF ( intermediate_timestep_count == 3 )  diss_prod3(:,j,i) = tend(:,j,i) - diss_adve3(:,j,i)
     2568
     2569       dum_pro = tend(:,j,i) - dum_adv                                          !> @todo remove later
     2570
     2571!
     2572!--    Diffusion term of TKE dissipation rate
    23202573       CALL diffusion_diss( i, j )
    23212574
     2575       IF ( intermediate_timestep_count == 1 )  diss_diff1(:,j,i) = tend(:,j,i) - diss_adve1(:,j,i) - diss_prod1(:,j,i) !> @todo remove later
     2576       IF ( intermediate_timestep_count == 2 )  diss_diff2(:,j,i) = tend(:,j,i) - diss_adve2(:,j,i) - diss_prod2(:,j,i)
     2577       IF ( intermediate_timestep_count == 3 )  diss_diff3(:,j,i) = tend(:,j,i) - diss_adve3(:,j,i) - diss_prod3(:,j,i)
     2578       IF ( intermediate_timestep_count == 3 )  dummy3(:,j,i) = km(:,j,i)
     2579
     2580       dum_dif = tend(:,j,i) - dum_adv - dum_pro                                !> @todo remove later
     2581
    23222582!
    23232583!--    Additional sink term for flows through plant canopies
    2324 !        IF ( plant_canopy )  CALL pcm_tendency( i, j, ? )                      !### not yet implemented
    2325 
    2326 !        CALL user_actions( i, j, 'diss-tendency' )                             !### not yet implemented
     2584!        IF ( plant_canopy )  CALL pcm_tendency( i, j, ? )                      !> @todo not yet implemented
     2585
     2586!        CALL user_actions( i, j, 'diss-tendency' )                             !> @todo not yet implemented
    23272587
    23282588!
     
    23382598                                                BTEST( wall_flags_0(k,j,i), 0 )&
    23392599                                               )
    2340           IF ( diss_p(k,j,i) <= 0.0_wp )  diss_p(k,j,i) = 0.1_wp * diss(k,j,i)
    23412600       ENDDO
    23422601
     
    23502609       ENDDO
    23512610
     2611       DO  l = 0, 1
     2612          surf_s = surf_def_v(l)%start_index(j,i)
     2613          surf_e = surf_def_v(l)%end_index(j,i)
     2614          DO  m = surf_s, surf_e
     2615             k = surf_def_v(l)%k(m)
     2616             diss_p(k,j,i) = surf_def_v(l)%us(m)**3 / ( kappa * 0.5_wp * dy )
     2617          ENDDO
     2618       ENDDO
     2619
     2620       DO  l = 2, 3
     2621          surf_s = surf_def_v(l)%start_index(j,i)
     2622          surf_e = surf_def_v(l)%end_index(j,i)
     2623          DO  m = surf_s, surf_e
     2624             k = surf_def_v(l)%k(m)
     2625             diss_p(k,j,i) = surf_def_v(l)%us(m)**3 / ( kappa * 0.5_wp * dx )
     2626          ENDDO
     2627       ENDDO
    23522628!
    23532629!--    Calculate tendencies for the next Runge-Kutta step
     
    23652641          ENDIF
    23662642       ENDIF
    2367 
    2368 !        IF ( intermediate_timestep_count == 1 )  dummy1(:,j,i) = e_p(:,j,i)
    2369 !        IF ( intermediate_timestep_count == 2 )  dummy2(:,j,i) = e_p(:,j,i)
    2370 !        IF ( intermediate_timestep_count == 3 )  dummy3(:,j,i) = e_p(:,j,i)
     2643!
     2644!--    Limit change of diss to be between -90% and +100%. Also, set an absolute
     2645!--    minimum value
     2646       DO  k = nzb, nzt+1
     2647          diss_p(k,j,i) = MIN( MAX( diss_p(k,j,i),         &
     2648                                    0.1_wp * diss(k,j,i),  &
     2649                                    0.0001_wp ),           &
     2650                               2.0_wp * diss(k,j,i) )
     2651       ENDDO
     2652
     2653       IF ( intermediate_timestep_count == 1 )  dummy1(:,j,i) = diss_p(:,j,i)   !> @todo remove later
     2654       IF ( intermediate_timestep_count == 2 )  dummy2(:,j,i) = diss_p(:,j,i)
     2655
     2656!        if ( i >= pis .and. i <= pie .and. j >= pjs .and. j <= pje ) then        !> @todo remove later
     2657!           WRITE(9, *) '---'
     2658!           WRITE(9, *) 'diss:'
     2659!           WRITE(9, 5555) 'adv   :', dum_adv(pkb:pkt)
     2660!           WRITE(9, 5555) 'pro   :', dum_pro(pkb:pkt)
     2661!           WRITE(9, 5555) 'dif   :', dum_dif(pkb:pkt)
     2662!           WRITE(9, 5555) 'tend  :', tend(pkb:pkt,j,i)
     2663!           WRITE(9, 5555) 'diss_p:', diss_p(pkb:pkt,j,i)
     2664!           WRITE(9, 5555) 'diss  :', diss(pkb:pkt,j,i)
     2665!           WRITE(9, *) '---'
     2666!           WRITE(9, 5555) 'km    :', km(pkb:pkt,j,i)
     2667!           flush(9)
     2668!        endif
    23712669
    23722670    ENDIF   ! dissipation equation
     
    23822680!> @warning The case with constant_flux_layer = F and use_surface_fluxes = T is
    23832681!>          not considered well!
     2682!> @todo Adjust production term in case of rans_tke_e simulation
    23842683!------------------------------------------------------------------------------!
    23852684 SUBROUTINE production_e
     
    30993398!> @warning The case with constant_flux_layer = F and use_surface_fluxes = T is
    31003399!>          not considered well!
     3400!> @todo non-neutral case is not yet considered for RANS mode
    31013401!------------------------------------------------------------------------------!
    3102  SUBROUTINE production_e_ij( i, j )
     3402 SUBROUTINE production_e_ij( i, j, diss_production )
    31033403
    31043404    USE arrays_3d,                                                             &
     
    31213421
    31223422    IMPLICIT NONE
     3423
     3424    LOGICAL :: diss_production
    31233425
    31243426    INTEGER(iwp) ::  i       !< running index x-direction
     
    31533455    REAL(wp), DIMENSION(nzb+1:nzt)  ::  dwdy        !< Gradient of w-component in y-direction
    31543456    REAL(wp), DIMENSION(nzb+1:nzt)  ::  dwdz        !< Gradient of w-component in z-direction
    3155 
     3457    REAL(wp), DIMENSION(nzb+1:nzt)  ::  tend_temp   !< temporal tendency
    31563458
    31573459    IF ( constant_flux_layer )  THEN
     
    33483650       ENDDO
    33493651
     3652!        IF ( .NOT. rans_tke_e )  THEN
     3653
    33503654       DO  k = nzb+1, nzt
    33513655
    3352           def = 2.0_wp * ( dudx(k)**2 + dvdy(k)**2 + dwdz(k)**2 ) +            &
    3353                            dudy(k)**2 + dvdx(k)**2 + dwdx(k)**2   +            &
    3354                            dwdy(k)**2 + dudz(k)**2 + dvdz(k)**2   +            &
    3355                 2.0_wp * ( dvdx(k)*dudy(k) + dwdx(k)*dudz(k) + dwdy(k)*dvdz(k) )
    3356 
    3357 !
    3358 !--       Production term according to Kato and Launder (1993)
    3359 !           def = SQRT( ( dudx(k) + dudy(k) + dudz(k) +                       &
    3360 !                         dvdx(k) + dvdy(k) + dvdz(k) +                       &
    3361 !                         dwdx(k) + dwdy(k) + dwdz(k)                         &
    3362 !                       )**4 -                                                &
    3363 !                       ( dudx(k)**2 + dvdy(k)**2 + dwdz(k)**2 +              &
    3364 !                         2.0_wp * ( dudy(k) * dvdx(k) +                      &
    3365 !                                    dudz(k) * dwdx(k) +                      &
    3366 !                                    dvdz(k) * dwdy(k) )                      &
    3367 !                       )**2                                                  &
    3368 !                     )
     3656          def = 2.0_wp * ( dudx(k)**2 + dvdy(k)**2 + dwdz(k)**2 ) +         &
     3657                           dudy(k)**2 + dvdx(k)**2 + dwdx(k)**2   +         &
     3658                           dwdy(k)**2 + dudz(k)**2 + dvdz(k)**2   +         &
     3659                2.0_wp * ( dvdx(k)*dudy(k) +                                &
     3660                           dwdx(k)*dudz(k) +                                &
     3661                           dwdy(k)*dvdz(k) )
    33693662
    33703663          IF ( def < 0.0_wp )  def = 0.0_wp
     
    33723665          flag  = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
    33733666
    3374           tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def * flag
     3667          tend_temp(k) = km(k,j,i) * def * flag
    33753668
    33763669       ENDDO
    33773670
     3671!        ELSE
     3672!
     3673!           DO  k = nzb+1, nzt
     3674! !
     3675! !--          Production term according to Kato and Launder (1993)
     3676!              def = SQRT( ( dudy(k)**2 + dvdz(k)**2 + dwdx(k)**2 +              &
     3677!                            dudz(k)**2 + dvdx(k)**2 + dwdy(k)**2 +              &
     3678!                            2.0_wp * ( dudy(k) * dvdx(k) +                      &
     3679!                                       dvdz(k) * dwdy(k) +                      &
     3680!                                       dwdx(k) * dudz(k) )       )              &
     3681!                        * ( dudy(k)**2 + dvdz(k)**2 + dwdx(k)**2 +              &
     3682!                            dudz(k)**2 + dvdx(k)**2 + dwdy(k)**2 -              &
     3683!                            2.0_wp * ( dudy(k) * dvdx(k) +                      &
     3684!                                       dvdz(k) * dwdy(k) +                      &
     3685!                                       dwdx(k) * dudz(k) )       )              &
     3686!                        )
     3687!
     3688!              IF ( def < 0.0_wp )  def = 0.0_wp
     3689!
     3690!              flag  = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
     3691!
     3692!              tend_temp(k) = km(k,j,i) * def * flag
     3693!
     3694!           ENDDO
     3695!
     3696!        ENDIF
     3697
     3698    ELSE  ! not constant_flux_layer
     3699
     3700!        IF ( .NOT. rans_tke_e )  THEN
     3701!
     3702!--       Calculate TKE production by shear. Here, no additional
     3703!--       wall-bounded code is considered.
     3704!--       Why?
     3705          DO  k = nzb+1, nzt
     3706
     3707             dudx(k)  =           ( u(k,j,i+1) - u(k,j,i)     ) * ddx
     3708             dudy(k)  = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) -                &
     3709                                    u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
     3710             dudz(k)  = 0.5_wp  * ( u(k+1,j,i) + u(k+1,j,i+1) -                &
     3711                                    u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
     3712
     3713             dvdx(k)  = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) -                &
     3714                                    v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
     3715             dvdy(k)  =           ( v(k,j+1,i) - v(k,j,i)     ) * ddy
     3716             dvdz(k)  = 0.5_wp  * ( v(k+1,j,i) + v(k+1,j+1,i) -                &
     3717                                    v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
     3718
     3719             dwdx(k)  = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) -                &
     3720                                    w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
     3721             dwdy(k)  = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) -                &
     3722                                    w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
     3723             dwdz(k)  =           ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
     3724
     3725             def = 2.0_wp * ( dudx(k)**2 + dvdy(k)**2 + dwdz(k)**2 ) +         &
     3726                              dudy(k)**2 + dvdx(k)**2 + dwdx(k)**2   +         &
     3727                              dwdy(k)**2 + dudz(k)**2 + dvdz(k)**2   +         &
     3728                   2.0_wp * ( dvdx(k)*dudy(k) +                                &
     3729                              dwdx(k)*dudz(k) +                                &
     3730                              dwdy(k)*dvdz(k) )
     3731
     3732             IF ( def < 0.0_wp )  def = 0.0_wp
     3733
     3734             flag        = MERGE( 1.0_wp, 0.0_wp,                              &
     3735                                  BTEST( wall_flags_0(k,j,i), 29 ) )
     3736             tend_temp(k) = km(k,j,i) * def * flag
     3737
     3738          ENDDO
     3739
     3740!        ELSE
     3741!
     3742!           DO  k = nzb+1, nzt
     3743!
     3744!              dudx(k)  =           ( u(k,j,i+1) - u(k,j,i)     ) * ddx
     3745!              dudy(k)  = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) -                &
     3746!                                     u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
     3747!              dudz(k)  = 0.5_wp  * ( u(k+1,j,i) + u(k+1,j,i+1) -                &
     3748!                                     u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
     3749!
     3750!              dvdx(k)  = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) -                &
     3751!                                     v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
     3752!              dvdy(k)  =           ( v(k,j+1,i) - v(k,j,i)     ) * ddy
     3753!              dvdz(k)  = 0.5_wp  * ( v(k+1,j,i) + v(k+1,j+1,i) -                &
     3754!                                     v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
     3755!
     3756!              dwdx(k)  = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) -                &
     3757!                                     w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
     3758!              dwdy(k)  = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) -                &
     3759!                                     w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
     3760!              dwdz(k)  =           ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
     3761! !
     3762! !--          Production term according to Kato and Launder (1993)
     3763!              def = SQRT( ( dudy(k)**2 + dvdz(k)**2 + dwdx(k)**2 +              &
     3764!                            dudz(k)**2 + dvdx(k)**2 + dwdy(k)**2 +              &
     3765!                            2.0_wp * ( dudy(k) * dvdx(k) +                      &
     3766!                                       dvdz(k) * dwdy(k) +                      &
     3767!                                       dwdx(k) * dudz(k) )       )              &
     3768!                        * ( dudy(k)**2 + dvdz(k)**2 + dwdx(k)**2 +              &
     3769!                            dudz(k)**2 + dvdx(k)**2 + dwdy(k)**2 -              &
     3770!                            2.0_wp * ( dudy(k) * dvdx(k) +                      &
     3771!                                       dvdz(k) * dwdy(k) +                      &
     3772!                                       dwdx(k) * dudz(k) )       )              &
     3773!                        )
     3774!
     3775!              IF ( def < 0.0_wp )  def = 0.0_wp
     3776!
     3777!              flag        = MERGE( 1.0_wp, 0.0_wp,                              &
     3778!                                   BTEST( wall_flags_0(k,j,i), 29 ) )
     3779!              tend_temp(k) = km(k,j,i) * def * flag
     3780!
     3781!           ENDDO
     3782!
     3783!        ENDIF
     3784
     3785    ENDIF
     3786
     3787    IF ( .NOT. diss_production )  THEN
     3788!
     3789!--    Production term in case of TKE production
     3790       DO  k = nzb+1, nzt
     3791          tend(k,j,i) = tend(k,j,i) + tend_temp(k)
     3792       ENDDO
    33783793    ELSE
    33793794!
    3380 !--    Calculate TKE production by shear. Here, no additional
    3381 !--    wall-bounded code is considered.
    3382 !--    Why?
     3795!--    Production term in case of dissipation-rate production (rans_tke_e)
    33833796       DO  k = nzb+1, nzt
    33843797
    3385           dudx(k)  =           ( u(k,j,i+1) - u(k,j,i)     ) * ddx
    3386           dudy(k)  = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) -                   &
    3387                                  u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
    3388           dudz(k)  = 0.5_wp  * ( u(k+1,j,i) + u(k+1,j,i+1) -                   &
    3389                                  u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
    3390 
    3391           dvdx(k)  = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) -                   &
    3392                                  v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
    3393           dvdy(k)  =           ( v(k,j+1,i) - v(k,j,i)     ) * ddy
    3394           dvdz(k)  = 0.5_wp  * ( v(k+1,j,i) + v(k+1,j+1,i) -                   &
    3395                                  v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
    3396 
    3397           dwdx(k)  = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) -                   &
    3398                                  w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
    3399           dwdy(k)  = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) -                   &
    3400                                  w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
    3401           dwdz(k)  =           ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
    3402 
    3403           def = 2.0_wp * ( dudx(k)**2 + dvdy(k)**2 + dwdz(k)**2 ) +            &
    3404                            dudy(k)**2 + dvdx(k)**2 + dwdx(k)**2   +            &
    3405                            dwdy(k)**2 + dudz(k)**2 + dvdz(k)**2   +            &
    3406                 2.0_wp * ( dvdx(k)*dudy(k) + dwdx(k)*dudz(k) + dwdy(k)*dvdz(k) )
    3407 
    3408 !
    3409 !--       Production term according to Kato and Launder (1993)
    3410 !           def = SQRT( ( dudx(k) + dudy(k) + dudz(k) +                        &
    3411 !                         dvdx(k) + dvdy(k) + dvdz(k) +                        &
    3412 !                         dwdx(k) + dwdy(k) + dwdz(k)                          &
    3413 !                       )**4 -                                                 &
    3414 !                       ( dudx(k)**2 + dvdy(k)**2 + dwdz(k)**2 +               &
    3415 !                         2.0_wp * ( dudy(k) * dvdx(k) +                       &
    3416 !                                    dudz(k) * dwdx(k) +                       &
    3417 !                                    dvdz(k) * dwdy(k) )                       &
    3418 !                       )**2                                                   &
    3419 !                     )
    3420 
    3421           IF ( def < 0.0_wp )  def = 0.0_wp
    3422 
    3423           flag        = MERGE( 1.0_wp, 0.0_wp,                                 &
    3424                                BTEST( wall_flags_0(k,j,i), 29 ) )
    3425           tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def * flag
    3426 
     3798          ! Standard TKE-e closure
     3799          tend(k,j,i) = tend(k,j,i) + tend_temp(k) * diss(k,j,i)               &
     3800                                    /( e(k,j,i) + 1.0E-20_wp )                 &
     3801                                    * c_1
     3802!           ! Production according to Koblitz (2013)
     3803!           tend(k,j,i) = tend(k,j,i) + tend_temp(k) * diss(k,j,i)               &
     3804!                                     /( e(k,j,i) + 1.0E-20_wp )                 &
     3805!                                     * ( c_1 + ( c_2 - c_1 )                    &
     3806!                                             * l_wall(k,j,i) / l_max )
     3807!           ! Production according to Detering and Etling (1985)
     3808!           !> @todo us is not correct if there are vertical walls
     3809!           tend(k,j,i) = tend(k,j,i) + tend_temp(k) * SQRT(e(k,j,i))            &
     3810!                                     * c_1 * c_0**3 / c_4 * f                   &
     3811!              / surf_def_h(0)%us(surf_def_h(0)%start_index(j,i))
    34273812       ENDDO
    3428 
    34293813    ENDIF
    34303814
     
    38224206    REAL(wp)     ::  l              !< mixing length
    38234207    REAL(wp)     ::  ll             !< adjusted l
    3824     REAL(wp)     ::  var_reference  !<
     4208    REAL(wp)     ::  var_reference  !< reference temperature
    38254209
    38264210#if defined( __nopointer )
    3827     REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  var  !<
     4211    REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  var  !< temperature
    38284212#else
    3829     REAL(wp), DIMENSION(:,:,:), POINTER ::  var  !<
     4213    REAL(wp), DIMENSION(:,:,:), POINTER ::  var  !< temperature
    38304214#endif
    38314215    REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) ::  dissipation  !< TKE dissipation
     
    38374221       DO  j = nys, nyn
    38384222          DO  k = nzb+1, nzt
     4223!
     4224!--          Predetermine flag to mask topography
     4225             flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
    38394226
    38404227!
     
    38514238                CALL mixing_length_rans( i, j, k, l, ll, var, var_reference )
    38524239
    3853                 dissipation(k,j) = c_m**3 * e(k,j,i) * SQRT( e(k,j,i) ) / ll
     4240                dissipation(k,j) = c_0**3 * e(k,j,i) * SQRT( e(k,j,i) ) / ll
     4241
     4242                diss(k,j,i) = dissipation(k,j) * flag
    38544243
    38554244             ELSEIF ( rans_tke_e )  THEN
     
    38594248             ENDIF
    38604249
    3861 !
    3862 !--          Predetermine flag to mask topography
    3863              flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
    3864 
    3865              tend(k,j,i) = tend(k,j,i)                                         &
    3866                                      + (                                       &
     4250             tend(k,j,i) = tend(k,j,i) + (                                     &
     4251                                           (                                   &
    38674252                       ( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) )     &
    38684253                     - ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) )     &
    3869                                        ) * ddx2  * flag                        &
    3870                                      + (                                       &
     4254                                           ) * ddx2  * flag                    &
     4255                                         + (                                   &
    38714256                       ( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) )     &
    38724257                     - ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) )     &
    3873                                        ) * ddy2  * flag                        &
    3874                                      + (                                       &
     4258                                           ) * ddy2  * flag                    &
     4259                                         + (                                   &
    38754260            ( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1)    &
    38764261                                                          * rho_air_zw(k)      &
    38774262          - ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k)      &
    38784263                                                          * rho_air_zw(k-1)    &
    3879                                        ) * ddzw(k) * drho_air(k) * flag        &
     4264                                           ) * ddzw(k) * drho_air(k)           &
     4265                                         ) * flag * dsig_e                     &
    38804266                          - dissipation(k,j) * flag
    38814267
     
    39594345    REAL(wp)     ::  l              !< mixing length
    39604346    REAL(wp)     ::  ll             !< adjusted l
    3961     REAL(wp)     ::  var_reference  !<
     4347    REAL(wp)     ::  var_reference  !< reference temperature
    39624348
    39634349#if defined( __nopointer )
    3964     REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  var  !<
     4350    REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  var  !< temperature
    39654351#else
    3966     REAL(wp), DIMENSION(:,:,:), POINTER ::  var     !<
     4352    REAL(wp), DIMENSION(:,:,:), POINTER ::  var     !< temperature
    39674353#endif
    39684354    REAL(wp), DIMENSION(nzb+1:nzt) ::  dissipation  !< dissipation of TKE
     
    39914377          CALL mixing_length_rans( i, j, k, l, ll, var, var_reference  )
    39924378
    3993           dissipation(k) = c_m**3 * e(k,j,i) * SQRT( e(k,j,i) ) / ll
     4379          dissipation(k) = c_0**3 * e(k,j,i) * SQRT( e(k,j,i) ) / ll
     4380
     4381          diss(k,j,i) = dissipation(k) * flag
    39944382
    39954383       ELSEIF ( rans_tke_e )  THEN
     
    40014389!
    40024390!--    Calculate the tendency term
    4003        tend(k,j,i) = tend(k,j,i)                                               &
    4004                                     + (                                        &
     4391       tend(k,j,i) = tend(k,j,i) + (                                           &
     4392                                      (                                        &
    40054393                      ( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) )      &
    40064394                    - ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) )      &
    4007                                       ) * ddx2  * flag / sig_e                 &
     4395                                      ) * ddx2                                 &
    40084396                                    + (                                        &
    40094397                      ( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) )      &
    40104398                    - ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) )      &
    4011                                       ) * ddy2  * flag / sig_e                 &
     4399                                      ) * ddy2                                 &
    40124400                                    + (                                        &
    40134401           ( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1)     &
     
    40154403         - ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k)       &
    40164404                                                         * rho_air_zw(k-1)     &
    4017                                       ) * ddzw(k) * drho_air(k) * flag / sig_e &
    4018                                     - dissipation(k) * flag
     4405                                      ) * ddzw(k) * drho_air(k)                &
     4406                                   ) * flag * dsig_e                           &
     4407                                 - dissipation(k) * flag
    40194408
    40204409    ENDDO
     
    40824471             flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
    40834472
    4084              tend(k,j,i) = tend(k,j,i)                                         &
    4085                                + (                                             &
     4473             tend(k,j,i) = tend(k,j,i) +                                       &
     4474                         (      (                                             &
    40864475                 ( km(k,j,i)+km(k,j,i+1) ) * ( diss(k,j,i+1)-diss(k,j,i) )     &
    40874476               - ( km(k,j,i)+km(k,j,i-1) ) * ( diss(k,j,i)-diss(k,j,i-1) )     &
    4088                                  ) * ddx2  * flag                              &
     4477                                 ) * ddx2                                      &
    40894478                               + (                                             &
    40904479                 ( km(k,j,i)+km(k,j+1,i) ) * ( diss(k,j+1,i)-diss(k,j,i) )     &
    40914480               - ( km(k,j,i)+km(k,j-1,i) ) * ( diss(k,j,i)-diss(k,j-1,i) )     &
    4092                                  ) * ddy2  * flag                              &
     4481                                 ) * ddy2                                      &
    40934482                               + (                                             &
    40944483      ( km(k,j,i)+km(k+1,j,i) ) * ( diss(k+1,j,i)-diss(k,j,i) ) * ddzu(k+1)    &
     
    40964485    - ( km(k,j,i)+km(k-1,j,i) ) * ( diss(k,j,i)-diss(k-1,j,i) ) * ddzu(k)      &
    40974486                                                    * rho_air_zw(k-1)          &
    4098                                  ) * ddzw(k) * drho_air(k) * flag              &
    4099                                - c_2 * diss(k,j,i)**2                          &
    4100                                      / ( e(k,j,i) + 1.0E-20_wp ) * flag
     4487                                 ) * ddzw(k) * drho_air(k)                     &
     4488                         ) * flag * dsig_diss                                  &
     4489                         - c_2 * diss(k,j,i)**2                                &
     4490                               / ( e(k,j,i) + 1.0E-20_wp ) * flag
    41014491
    41024492          ENDDO
     
    41234513    IMPLICIT NONE
    41244514
    4125     INTEGER(iwp) ::  i              !< running index x direction
    4126     INTEGER(iwp) ::  j              !< running index y direction
    4127     INTEGER(iwp) ::  k              !< running index z direction
    4128 
    4129     REAL(wp)     ::  flag           !< flag to mask topography
    4130 
    4131     REAL(wp), DIMENSION(nzb+1:nzt) :: tend_temp
     4515    INTEGER(iwp) ::  i         !< running index x direction
     4516    INTEGER(iwp) ::  j         !< running index y direction
     4517    INTEGER(iwp) ::  k         !< running index z direction
     4518
     4519    REAL(wp)     ::  flag      !< flag to mask topography
    41324520
    41334521!
     
    41414529!
    41424530!--    Calculate the tendency term
    4143        tend_temp(k) =           (                                              &
     4531       tend(k,j,i) =  tend(k,j,i) +                                            &
     4532                   (            (                                              &
    41444533                ( km(k,j,i)+km(k,j,i+1) ) * ( diss(k,j,i+1)-diss(k,j,i) )      &
    41454534              - ( km(k,j,i)+km(k,j,i-1) ) * ( diss(k,j,i)-diss(k,j,i-1) )      &
    4146                                 ) * ddx2  * flag / sig_diss                    &
     4535                                ) * ddx2                                       &
    41474536                              + (                                              &
    41484537                ( km(k,j,i)+km(k,j+1,i) ) * ( diss(k,j+1,i)-diss(k,j,i) )      &
    41494538              - ( km(k,j,i)+km(k,j-1,i) ) * ( diss(k,j,i)-diss(k,j-1,i) )      &
    4150                                 ) * ddy2  * flag / sig_diss                    &
     4539                                ) * ddy2                                       &
    41514540                              + (                                              &
    41524541     ( km(k,j,i)+km(k+1,j,i) ) * ( diss(k+1,j,i)-diss(k,j,i) ) * ddzu(k+1)     &
     
    41544543   - ( km(k,j,i)+km(k-1,j,i) ) * ( diss(k,j,i)-diss(k-1,j,i) ) * ddzu(k)       &
    41554544                                                   * rho_air_zw(k-1)           &
    4156                                 ) * ddzw(k) * drho_air(k) * flag / sig_diss    &
    4157                               - c_2 * diss(k,j,i)**2                           &
    4158                                     / ( e(k,j,i) + 1.0E-20_wp ) * flag
    4159 
    4160        tend(k,j,i) = tend(k,j,i) + tend_temp(k)
     4545                                ) * ddzw(k) * drho_air(k)                      &
     4546                   ) * flag * dsig_diss                                        &
     4547                   - c_2 * diss(k,j,i)**2 / ( e(k,j,i) + 1.0E-20_wp ) * flag
    41614548
    41624549    ENDDO
     
    42904677!> Computation of the turbulent diffusion coefficients for momentum and heat
    42914678!> according to Prandtl-Kolmogorov.
     4679!> @todo consider non-default surfaces
    42924680!------------------------------------------------------------------------------!
    42934681 SUBROUTINE tcm_diffusivities( var, var_reference )
     
    42974685        ONLY:  e_min, outflow_l, outflow_n, outflow_r, outflow_s
    42984686
     4687    USE grid_variables,                                                        &
     4688        ONLY:  dx, dy
     4689
    42994690    USE statistics,                                                            &
    43004691        ONLY :  rmask, sums_l_l
    43014692
    43024693    USE surface_mod,                                                           &
    4303         ONLY :  bc_h, surf_def_h
     4694        ONLY :  bc_h, surf_def_h, surf_def_v
    43044695
    43054696    IMPLICIT NONE
    43064697
    4307     INTEGER(iwp) ::  i                   !<
    4308     INTEGER(iwp) ::  j                   !<
    4309     INTEGER(iwp) ::  k                   !<
    4310     INTEGER(iwp) ::  m                   !<
    4311     INTEGER(iwp) ::  n                   !<
    4312     INTEGER(iwp) ::  omp_get_thread_num  !<
    4313     INTEGER(iwp) ::  sr                  !<
    4314     INTEGER(iwp) ::  tn                  !<
    4315 
    4316     REAL(wp)     ::  flag                !<
    4317     REAL(wp)     ::  l                   !<
    4318     REAL(wp)     ::  ll                  !<
    4319     REAL(wp)     ::  var_reference       !<
     4698    INTEGER(iwp) ::  i                   !< loop index
     4699    INTEGER(iwp) ::  j                   !< loop index
     4700    INTEGER(iwp) ::  k                   !< loop index
     4701    INTEGER(iwp) ::  m                   !< loop index
     4702    INTEGER(iwp) ::  n                   !< loop index
     4703    INTEGER(iwp) ::  omp_get_thread_num  !< opemmp function to get thread number
     4704    INTEGER(iwp) ::  sr                  !< statistic region
     4705    INTEGER(iwp) ::  tn                  !< thread number
     4706
     4707    REAL(wp)     ::  flag                !< topography flag
     4708    REAL(wp)     ::  l                   !< mixing length
     4709    REAL(wp)     ::  ll                  !< adjusted mixing length
     4710    REAL(wp)     ::  var_reference       !< reference temperature
    43204711
    43214712#if defined( __nopointer )
    4322     REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  var  !<
     4713    REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  var  !< temperature
    43234714#else
    4324     REAL(wp), DIMENSION(:,:,:), POINTER ::  var  !<
     4715    REAL(wp), DIMENSION(:,:,:), POINTER ::  var  !< temperature
    43254716#endif
    43264717
     
    43654756!
    43664757!--             Compute diffusion coefficients for momentum and heat
    4367                 km(k,j,i) = c_m * l * SQRT( e(k,j,i) ) * flag
     4758                km(k,j,i) = c_0 * l * SQRT( e(k,j,i) ) * flag
    43684759                kh(k,j,i) = ( 1.0_wp + 2.0_wp * l / ll ) * km(k,j,i) * flag
    43694760!
     
    43914782!
    43924783!--             Compute diffusion coefficients for momentum and heat
    4393                 km(k,j,i) = c_m * l * SQRT( e(k,j,i) ) * flag
     4784                km(k,j,i) = c_0 * l * SQRT( e(k,j,i) ) * flag
    43944785                kh(k,j,i) = km(k,j,i) / prandtl_number * flag
    43954786!
     
    44144805!
    44154806!--             Compute diffusion coefficients for momentum and heat
    4416                 km(k,j,i) = c_mu * e(k,j,i)**2 / ( diss(k,j,i) + 1.E-10 ) * flag
     4807                km(k,j,i) = c_0**4 * e(k,j,i)**2 / ( diss(k,j,i) + 1.0E-30_wp ) * flag
    44174808                kh(k,j,i) = km(k,j,i) / prandtl_number * flag
     4809!
     4810!--             Summation for averaged profile of mixing length (cf. flow_statistics)
     4811                DO  sr = 0, statistic_regions
     4812                   sums_l_l(k,sr,tn) = sums_l_l(k,sr,tn) +                     &
     4813                      c_0**3 * e(k,j,i) * SQRT(e(k,j,i)) /                     &
     4814                      ( diss(k,j,i) + 1.0E-30_wp ) * rmask(j,i,sr) * flag
     4815                ENDDO
    44184816
    44194817             ENDDO
     
    44344832!-- so far vertical surfaces require usage of a Prandtl-layer where the boundary
    44354833!-- values of the diffusivities are not needed.
    4436 
    44374834    IF ( .NOT. rans_tke_e )  THEN
    44384835!
     
    44684865          ENDDO
    44694866       ENDDO
     4867!
     4868!--    North- and southward facing surfaces
     4869       DO  n = 0, 1
     4870          DO  m = 1, surf_def_v(n)%ns
     4871             i = surf_def_v(n)%i(m)
     4872             j = surf_def_v(n)%j(m)
     4873             k = surf_def_v(n)%k(m)
     4874             km(k,j,i) = kappa * surf_def_v(n)%us(m) * 0.5_wp * dy
     4875             kh(k,j,i) = 1.35_wp * km(k,j,i)
     4876          ENDDO
     4877       ENDDO
     4878!
     4879!--    West- and eastward facing surfaces
     4880       DO  n = 2, 3
     4881          DO  m = 1, surf_def_v(n)%ns
     4882             i = surf_def_v(n)%i(m)
     4883             j = surf_def_v(n)%j(m)
     4884             k = surf_def_v(n)%k(m)
     4885             km(k,j,i) = kappa * surf_def_v(n)%us(m) * 0.5_wp * dx
     4886             kh(k,j,i) = 1.35_wp * km(k,j,i)
     4887          ENDDO
     4888       ENDDO
    44704889
    44714890       CALL exchange_horiz( km, nbgp )
     
    44734892
    44744893    ENDIF
     4894
    44754895!
    44764896!-- Model top
     
    45184938    INTEGER(iwp) ::  j      !< loop index y direction
    45194939    INTEGER(iwp) ::  k      !< loop index z direction
    4520     INTEGER, INTENT(IN) ::  mod_count  !<
     4940    INTEGER, INTENT(IN) ::  mod_count  !< flag defining where pointers point to
    45214941
    45224942#if defined( __nopointer )
Note: See TracChangeset for help on using the changeset viewer.