Changeset 3083 for palm/trunk/SOURCE


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

merge with branch rans

Location:
palm/trunk
Files:
14 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk

  • palm/trunk/SOURCE

  • palm/trunk/SOURCE/Makefile

    r2963 r3083  
    2525# -----------------
    2626# $Id$
     27# add turbulence_closure_mod to parin
     28#
     29# 2963 2018-04-12 14:47:44Z suehring
    2730# Introduce index for vegetation/wall, pavement/green-wall and water/window
    2831# surfaces, for clearer access of surface fraction, albedo, emissivity, etc. .
     
    12301233        spectra_mod.o \
    12311234        synthetic_turbulence_generator_mod.o \
     1235        turbulence_closure_mod.o \
    12321236        uv_exposure_model_mod.o \
    12331237        vertical_nesting_mod.o \
  • palm/trunk/SOURCE/check_parameters.f90

    r3065 r3083  
    2424! Former revisions:
    2525! -----------------
    26 ! $Id$
     26! $Id: check_parameters.f90 2520 2017-10-05 13:50:26Z gronemeier &
     27! Add inital profile output for e (TG)
     28!
     29! 3065 2018-06-12 07:03:02Z Giersch
    2730! dz was replaced by dz(1), error message revised
    2831!
     
    24552458             ENDIF
    24562459
    2457           CASE ( 'e' )
     2460          CASE ( 'e', '#e' )
    24582461             dopr_index(i)  = 8
    24592462             dopr_unit(i)   = 'm2/s2'
    24602463             hom(:,2,8,:)   = SPREAD( zu, 2, statistic_regions+1 )
    24612464             hom(nzb,2,8,:) = 0.0_wp
     2465             IF ( data_output_pr(i)(1:1) == '#' )  THEN
     2466                dopr_initial_index(i) = 8
     2467                hom(:,2,8,:)          = SPREAD( zu, 2, statistic_regions+1 )
     2468                data_output_pr(i)     = data_output_pr(i)(2:)
     2469             ENDIF
    24622470
    24632471          CASE ( 'km', '#km' )
  • palm/trunk/SOURCE/data_output_mask.f90

    r3045 r3083  
    2525! -----------------
    2626! $Id$
     27!
     28!
     29! 3045 2018-05-28 07:55:41Z Giersch
    2730! Error messages revised
    2831!
  • palm/trunk/SOURCE/header.f90

    r3065 r3083  
    2525! -----------------
    2626! $Id$
     27! Print RANS-mode constants
     28!
     29! 3065 2018-06-12 07:03:02Z Giersch
    2730! Header output concerning stretching revised
    2831!
     
    450453    USE synthetic_turbulence_generator_mod,                                    &
    451454        ONLY:  stg_header
     455
     456    USE turbulence_closure_mod,                                                &
     457        ONLY:  rans_const_c, rans_const_sigma
    452458
    453459    IMPLICIT NONE
     
    19461952       IF ( wall_adjustment )  WRITE ( io, 453 )  wall_adjustment_factor
    19471953    ENDIF
    1948 
     1954    IF ( rans_mode )  THEN
     1955       WRITE ( io, 457 )  rans_const_c, rans_const_sigma
     1956    ENDIF
    19491957!
    19501958!-- Special actions during the run
     
    23882396456 FORMAT ('    Day of the year at model start :   day_init = ',I3             &
    23892397            /'    UTC time at model start        :   time_utc_init = ',F7.1' s')
     2398457 FORMAT ('    RANS-mode constants: c_0 = ',F9.5/         &
     2399            '                         c_1 = ',F9.5/         &
     2400            '                         c_2 = ',F9.5/         &
     2401            '                         c_3 = ',F9.5/         &
     2402            '                         c_4 = ',F9.5/         &
     2403            '                         sigma_e    = ',F9.5/  &
     2404            '                         sigma_diss = ',F9.5)
    23902405470 FORMAT (//' Actions during the simulation:'/ &
    23912406              ' -----------------------------'/)
  • palm/trunk/SOURCE/init_3d_model.f90

  • palm/trunk/SOURCE/model_1d_mod.f90

    r3049 r3083  
    2525! -----------------
    2626! $Id$
     27! Bugfixes:
     28!   - preset te_diss and te_e to avoid runtime errors
     29!   - implementation of buoyancy term to dissipation
     30!     according to Sogachev et al. (2012)
     31!   - where diss_p < 0 set diss_p = 0.1 diss
     32!   - calculate progn eq(diss) starting from nzb+1
     33! Changes:
     34!   - add sig_e to TKE equation
     35!   - adjust prognostic equation of diss
     36!   - set model constants according to Koblitz (2013)
     37!   - renamed c_m to c_0
     38!   - rename l_black into l1d_init
     39!   - calculate l_grid within init_1d_model and save it as l1d_init
     40!   - calculate l1d according to DE85 if dissipation is a prognostic value
     41!   - made annotations doxygen-readable
     42!
     43! 3049 2018-05-29 13:52:36Z Giersch
    2744! Error messages revised
    2845!
     
    164181    USE kinds
    165182
    166     USE pegrid
     183    USE pegrid,                                                                &
     184        ONLY:  myid
    167185       
    168186
     
    175193    LOGICAL ::  stop_dt_1d = .FALSE.             !< termination flag, used in case of too small timestep (1d-model)
    176194
    177     REAL(wp) ::  c_1 = 1.44_wp                 !< model constant
    178     REAL(wp) ::  c_2 = 1.92_wp                 !< model constant
    179     REAL(wp) ::  c_3 = 1.44_wp                 !< model constant
    180     REAL(wp) ::  c_h = 0.0015_wp               !< model constant according to Detering and Etling (1985)
    181     REAL(wp) ::  c_m = 0.4_wp                  !< model constant, 0.4 according to Detering and Etling (1985)
    182     REAL(wp) ::  c_mu = 0.09_wp                !< model constant
     195    REAL(wp) ::  alpha_buoyancy                !< model constant according to Koblitz (2013)
     196    REAL(wp) ::  c_0 = 0.03_wp**0.25_wp        !< model constant according to Koblitz (2013)
     197    REAL(wp) ::  c_1 = 1.52_wp                 !< model constant according to Koblitz (2013)
     198    REAL(wp) ::  c_2 = 1.83_wp                 !< model constant according to Koblitz (2013)
     199    REAL(wp) ::  c_3                           !< model constant
     200    REAL(wp) ::  c_mu                          !< model constant
    183201    REAL(wp) ::  damp_level_1d = -1.0_wp       !< namelist parameter
    184202    REAL(wp) ::  dt_1d = 60.0_wp               !< dynamic timestep (1d-model)
     
    187205    REAL(wp) ::  dt_run_control_1d = 60.0_wp   !< namelist parameter
    188206    REAL(wp) ::  end_time_1d = 864000.0_wp     !< namelist parameter
     207    REAL(wp) ::  lambda                        !< maximum mixing length
    189208    REAL(wp) ::  qs1d                          !< characteristic humidity scale (1d-model)
    190209    REAL(wp) ::  simulated_time_1d = 0.0_wp    !< updated simulated time (1d-model)
    191     REAL(wp) ::  sig_diss = 1.3_wp             !< model constant
     210    REAL(wp) ::  sig_diss = 2.95_wp            !< model constant according to Koblitz (2013)
     211    REAL(wp) ::  sig_e = 2.95_wp               !< model constant according to Koblitz (2013)
    192212    REAL(wp) ::  time_pr_1d = 0.0_wp           !< updated simulated time for profile output (1d-model)
    193213    REAL(wp) ::  time_run_control_1d = 0.0_wp  !< updated simulated time for run-control output (1d-model)
     
    198218    REAL(wp) ::  z01d                          !< roughness length for momentum (1d-model)
    199219    REAL(wp) ::  z0h1d                         !< roughness length for scalars (1d-model)
    200 
    201220
    202221    REAL(wp), DIMENSION(:), ALLOCATABLE ::  diss1d   !< tke dissipation rate (1d-model)
     
    279298
    280299    INTEGER(iwp) ::  k  !< loop index
    281    
    282     REAL(wp) ::  lambda !< maximum mixing length
    283300
    284301!
     
    324341!
    325342!--       Use the same mixing length as in 3D model (LES-mode)
    326           !@todo: rename (delete?) this option
    327           ! As the mixing length is different between RANS and LES mode, it
    328           ! must be distinguished here between these modes. For RANS mode,
    329           ! the mixing length is calculated accoding to Blackadar, which is
    330           ! the other option at this point.
    331           ! Maybe delete this option entirely (not appropriate in LES case)
    332           ! 2018-03-20, gronemeier
     343          !> @todo rename (delete?) this option
     344          !> As the mixing length is different between RANS and LES mode, it
     345          !> must be distinguished here between these modes. For RANS mode,
     346          !> the mixing length is calculated accoding to Blackadar, which is
     347          !> the other option at this point.
     348          !> Maybe delete this option entirely (not appropriate in LES case)
     349          !> 2018-03-20, gronemeier
    333350          DO  k = nzb+1, nzt
    334351             l1d_init(k)  = ( dx * dy * dzw(k) )**0.33333333333333_wp
     
    359376       us1d = 0.1_wp   ! without initial friction the flow would not change
    360377    ELSE
    361        diss1d(nzb+1) = 1.0_wp
     378       diss1d(nzb+1) = 0.001_wp
    362379       e1d(nzb+1)  = 1.0_wp
    363380       km1d(nzb+1) = 1.0_wp
     
    372389
    373390!
    374 !-- Tendencies must be preset in order to avoid runtime errors within the
    375 !-- first Runge-Kutta step
     391!-- Tendencies must be preset in order to avoid runtime errors
     392    te_diss  = 0.0_wp
    376393    te_dissm = 0.0_wp
     394    te_e  = 0.0_wp
    377395    te_em = 0.0_wp
    378396    te_um = 0.0_wp
     
    381399!
    382400!-- Set model constant
    383     IF ( dissipation_1d == 'as_in_3d_model' )  c_m = 0.1_wp
     401    IF ( dissipation_1d == 'as_in_3d_model' )  c_0 = 0.1_wp
     402    c_mu = c_0**4
    384403
    385404!
     
    423442!-- Determine the time step at the start of a 1D-simulation and
    424443!-- determine and printout quantities used for run control
    425     dt_1d = 1.0_wp
     444    dt_1d = 0.01_wp
    426445    CALL run_control_1d
    427446
     
    483502!--             dissipation rate
    484503                IF ( dissipation_1d == 'detering' )  THEN
    485                    diss1d(k) = c_m**3 * e1d(k) * SQRT( e1d(k) ) / l1d_diss(k)
     504                   diss1d(k) = c_0**3 * e1d(k) * SQRT( e1d(k) ) / l1d_diss(k)
    486505                ELSEIF ( dissipation_1d == 'as_in_3d_model' )  THEN
    487506                   diss1d(k) = ( 0.19_wp + 0.74_wp * l1d_diss(k) / l1d_init(k) &
     
    497516                                     kmzp * ( e1d(k+1) - e1d(k) ) * ddzu(k+1)  &
    498517                                   - kmzm * ( e1d(k) - e1d(k-1) ) * ddzu(k)    &
    499                                                  ) * ddzw(k)                   &
     518                                                 ) * ddzw(k) / sig_e           &
    500519                                   - diss1d(k)
    501520
     
    503522!
    504523!--                dissipation rate
    505                    te_diss(k) = km1d(k) *                                      &
     524                   IF ( rif1d(k) >= 0.0_wp )  THEN
     525                      alpha_buoyancy = 1.0_wp - l1d(k) / lambda
     526                   ELSE
     527                      alpha_buoyancy = 1.0_wp - ( 1.0_wp + ( c_2 - 1.0_wp )    &
     528                                                         / ( c_2 - c_1    ) )  &
     529                                              * l1d(k) / lambda
     530                   ENDIF
     531                   c_3 = ( c_1 - c_2 ) * alpha_buoyancy
     532                   te_diss(k) = ( km1d(k) *                                    &
    506533                                  ( ( ( u1d(k+1) - u1d(k-1) ) * dd2zu(k) )**2  &
    507534                                  + ( ( v1d(k+1) - v1d(k-1) ) * dd2zu(k) )**2  &
    508                                   ) * c_1 * c_mu**0.75 / c_h * f / us1d        &
    509                                     * SQRT(e1d(k))                             &
     535                                  ) * ( c_1 + (c_2 - c_1) * l1d(k) / lambda )  &
    510536                                  - g / pt_0 * kh1d(k) * flux * c_3            &
    511                                     * diss1d(k) / ( e1d(k) + 1.0E-20_wp )      &
    512                                   + ( kmzp * ( diss1d(k+1) - diss1d(k) )       &
     537                                  - c_2 * diss1d(k)                            &
     538                                ) * diss1d(k) / ( e1d(k) + 1.0E-20_wp )        &
     539                                + (   kmzp * ( diss1d(k+1) - diss1d(k) )       &
    513540                                           * ddzu(k+1)                         &
    514541                                    - kmzm * ( diss1d(k) - diss1d(k-1) )       &
    515542                                           * ddzu(k)                           &
    516                                     ) * ddzw(k) / sig_diss                     &
    517                                   - c_2 * diss1d(k)**2 / ( e1d(k) + 1.0E-20_wp )
     543                                  ) * ddzw(k) / sig_diss
    518544
    519545                ENDIF
     
    546572!--          dissipation rate
    547573             IF ( dissipation_1d == 'detering' )  THEN
    548                 diss1d(k) = c_m**3 * e1d(k) * SQRT( e1d(k) ) / l1d_diss(k)
     574                diss1d(k) = c_0**3 * e1d(k) * SQRT( e1d(k) ) / l1d_diss(k)
    549575             ELSEIF ( dissipation_1d == 'as_in_3d_model' )  THEN
    550576                diss1d(k) = ( 0.19_wp + 0.74_wp * l1d_diss(k) / l1d_init(k) )  &
     
    565591!--          TKE
    566592             IF ( .NOT. dissipation_1d == 'prognostic' )  THEN
     593                !> @query why integrate over 2dz
     594                !>   Why is it allowed to integrate over two delta-z for e
     595                !>   while for u and v it is not?
     596                !>   2018-04-23, gronemeier
    567597                te_e(k) = km1d(k) * ( ( ( u1d(k+1) - u1d(k-1) ) * dd2zu(k) )**2&
    568598                                    + ( ( v1d(k+1) - v1d(k-1) ) * dd2zu(k) )**2&
     
    572602                                     kmzp * ( e1d(k+1) - e1d(k) ) * ddzu(k+1)  &
    573603                                   - kmzm * ( e1d(k) - e1d(k-1) ) * ddzu(k)    &
    574                                                  ) * ddzw(k)                   &
     604                                                 ) * ddzw(k) / sig_e           &
    575605                                   - diss1d(k)
    576606             ENDIF
     
    595625             ENDDO
    596626
    597              IF ( dissipation_1d == 'prognostic' )  THEN
    598                 DO  k = nzb_diff, nzt
    599                    diss1d_p(k) = diss1d(k) + dt_1d * ( tsc(2) * te_diss(k) + &
    600                                                  tsc(3) * te_dissm(k) )
    601                 ENDDO
    602              ENDIF
    603627!
    604628!--          Eliminate negative TKE values, which can result from the
     
    606630!--          value is reduced to 10 percent of its old value.
    607631             WHERE ( e1d_p < 0.0_wp )  e1d_p = 0.1_wp * e1d
     632
     633             IF ( dissipation_1d == 'prognostic' )  THEN
     634                DO  k = nzb+1, nzt
     635                   diss1d_p(k) = diss1d(k) + dt_1d * ( tsc(2) * te_diss(k) + &
     636                                                       tsc(3) * te_dissm(k) )
     637                ENDDO
     638                WHERE ( diss1d_p < 0.0_wp )  diss1d_p = 0.1_wp * diss1d
     639             ENDIF
    608640          ENDIF
    609641
     
    643675                   IF ( dissipation_1d == 'prognostic' )  THEN
    644676                      DO k = nzb+1, nzt
    645                          te_dissm(k) = -9.5625_wp * te_diss(k) + 5.3125_wp * te_dissm(k)
     677                         te_dissm(k) = -9.5625_wp * te_diss(k)  &
     678                                     +  5.3125_wp * te_dissm(k)
    646679                      ENDDO
    647680                   ENDIF
     
    650683             ENDIF
    651684          ENDIF
    652 
    653685
    654686!
     
    701733
    702734             ENDIF    ! constant_flux_layer
    703 
     735             !> @todo combine if clauses
     736             !>   The previous and following if clauses can be combined into a
     737             !>   single clause
     738             !>   2018-04-23, gronemeier
    704739!
    705740!--          Compute the Richardson-flux numbers,
     
    789824!--             compatibility with the 3D model.
    790825                IF ( ibc_e_b == 2 )  THEN
    791                    e1d(nzb+1) = ( us1d / c_m )**2
     826                   e1d(nzb+1) = ( us1d / c_0 )**2
    792827                ENDIF
    793828                IF ( dissipation_1d == 'prognostic' )  THEN
    794                    e1d(nzb+1) = us1d**2 / SQRT( c_mu )
     829                   e1d(nzb+1) = ( us1d / c_0 )**2
    795830                   diss1d(nzb+1) = us1d**3 / ( kappa * zu(nzb+1) )
    796831                   diss1d(nzb) = diss1d(nzb+1)
     
    829864!--          in the dissipation of TKE via l1d_diss. Otherwise, km1d would be
    830865!--          too large.
    831              IF ( mixing_length_1d == 'blackadar' )  THEN
     866             IF ( dissipation_1d /= 'prognostic' )  THEN
     867                IF ( mixing_length_1d == 'blackadar' )  THEN
     868                   DO  k = nzb+1, nzt
     869                      IF ( rif1d(k) >= 0.0_wp )  THEN
     870                         l1d(k) = l1d_init(k) / ( 1.0_wp + 5.0_wp * rif1d(k) )
     871                         l1d_diss(k) = l1d(k)
     872                      ELSE
     873                         l1d(k) = l1d_init(k)
     874                         l1d_diss(k) = l1d_init(k) *                           &
     875                                       SQRT( 1.0_wp - 16.0_wp * rif1d(k) )
     876                      ENDIF
     877                   ENDDO
     878                ELSEIF ( mixing_length_1d == 'as_in_3d_model' )  THEN
     879                   DO  k = nzb+1, nzt
     880                      dpt_dz = ( pt_init(k+1) - pt_init(k-1) ) * dd2zu(k)
     881                      IF ( dpt_dz > 0.0_wp )  THEN
     882                         l_stable = 0.76_wp * SQRT( e1d(k) )                   &
     883                                  / SQRT( g / pt_init(k) * dpt_dz ) + 1E-5_wp
     884                      ELSE
     885                         l_stable = l1d_init(k)
     886                      ENDIF
     887                      l1d(k) = MIN( l1d_init(k), l_stable )
     888                      l1d_diss(k) = l1d(k)
     889                   ENDDO
     890                ENDIF
     891             ELSE
    832892                DO  k = nzb+1, nzt
    833                    IF ( rif1d(k) >= 0.0_wp )  THEN
    834                       l1d(k) = l1d_init(k) / ( 1.0_wp + 5.0_wp * rif1d(k) )
    835                       l1d_diss(k) = l1d(k)
    836                    ELSE
    837                       l1d(k) = l1d_init(k)
    838                       l1d_diss(k) = l1d_init(k) *                              &
    839                                     SQRT( 1.0_wp - 16.0_wp * rif1d(k) )
    840                    ENDIF
    841                 ENDDO
    842              ELSEIF ( mixing_length_1d == 'as_in_3d_model' )  THEN
    843                 DO  k = nzb+1, nzt
    844                    dpt_dz = ( pt_init(k+1) - pt_init(k-1) ) * dd2zu(k)
    845                    IF ( dpt_dz > 0.0_wp )  THEN
    846                       l_stable = 0.76_wp * SQRT( e1d(k) ) /                    &
    847                                      SQRT( g / pt_init(k) * dpt_dz ) + 1E-5_wp
    848                    ELSE
    849                       l_stable = l1d_init(k)
    850                    ENDIF
    851                    l1d(k) = MIN( l1d_init(k), l_stable )
    852                    l1d_diss(k) = l1d(k)
     893                   l1d(k) = c_0**3 * e1d(k) * SQRT( e1d(k) )                   &
     894                          / ( diss1d(k) + 1.0E-30_wp )
    853895                ENDDO
    854896             ENDIF
     
    867909                ENDIF
    868910             ENDIF
     911
    869912             IF ( dissipation_1d == 'prognostic' )  THEN
    870913                DO  k = nzb_diff, nzt
     
    873916             ELSE
    874917                DO  k = nzb_diff, nzt
    875                    km1d(k) = c_m * SQRT( e1d(k) ) * l1d(k)
     918                   km1d(k) = c_0 * SQRT( e1d(k) ) * l1d(k)
    876919                ENDDO
    877920             ENDIF
     
    10181061
    10191062    REAL(wp) ::  dt_diff  !< time step accorind to diffusion criterion
     1063    REAL(wp) ::  dt_old   !< previous time step
    10201064    REAL(wp) ::  fac      !< factor of criterion
    10211065    REAL(wp) ::  value    !< auxiliary variable
    10221066
    10231067!
     1068!-- Save previous time step
     1069    dt_old = dt_1d
     1070
     1071!
    10241072!-- Compute the currently feasible time step according to the diffusion
    10251073!-- criterion. At nzb+1 the half grid length is used.
    1026     fac = 0.125  !0.35_wp                                                       !### changed from 0.35
     1074    fac = 0.125
    10271075    dt_diff = dt_max_1d
    10281076    DO  k = nzb+2, nzt
     
    10341082
    10351083!
     1084!-- Limit the new time step to a maximum of 10 times the previous time step
     1085    dt_1d = MIN( dt_old * 10.0_wp, dt_1d )
     1086
     1087!
    10361088!-- Set flag when the time step becomes too small
    1037     IF ( dt_1d < ( 0.00001_wp * dt_max_1d ) )  THEN
     1089    IF ( dt_1d < ( 1.0E-15_wp * dt_max_1d ) )  THEN
    10381090       stop_dt_1d = .TRUE.
    10391091
  • palm/trunk/SOURCE/modules.f90

    r3065 r3083  
    2525! -----------------
    2626! $Id$
     27! set dt_3d = 0.01
     28!
     29! 3065 2018-06-12 07:03:02Z Giersch
    2730! Variables concerning stretching introduced or revised
    2831!
     
    14221425    REAL(wp) ::  dt_run_control = 60.0_wp                      !< namelist parameter
    14231426    REAL(wp) ::  dt_spinup = 60.0_wp                           !< namelist parameter
    1424     REAL(wp) ::  dt_3d = 1.0_wp                                !< time step
     1427    REAL(wp) ::  dt_3d = 0.01_wp                               !< time step
    14251428    REAL(wp) ::  dz_max = 1000.0_wp                            !< namelist parameter
    14261429    REAL(wp) ::  dz_stretch_factor = 1.08_wp                   !< namelist parameter
  • palm/trunk/SOURCE/parin.f90

    r3065 r3083  
    2525! -----------------
    2626! $Id$
     27! Added rans_const_c and rans_const_sigma as input parameters (TG)
     28!
     29! 3065 2018-06-12 07:03:02Z Giersch
    2730! New initialization parameters added
    2831!
     
    467470    USE synthetic_turbulence_generator_mod,                                    &
    468471        ONLY:  stg_parin
     472
     473    USE turbulence_closure_mod,                                                &
     474        ONLY:  rans_const_c, rans_const_sigma
    469475
    470476    USE urban_surface_mod,                                                     &
     
    531537             pt_vertical_gradient_level, q_surface, q_surface_initial_change,  &
    532538             q_vertical_gradient, q_vertical_gradient_level,                   &
    533              random_generator, random_heatflux, rans_mode,                     &
     539             random_generator, random_heatflux, rans_const_c, rans_const_sigma,&
     540             rans_mode,                                                        &
    534541             rayleigh_damping_factor, rayleigh_damping_height,                 &
    535542             recycling_width, recycling_yshift,                                &
     
    601608             pt_vertical_gradient_level, q_surface, q_surface_initial_change,  &
    602609             q_vertical_gradient, q_vertical_gradient_level,                   &
    603              random_generator, random_heatflux, rans_mode,                     &
     610             random_generator, random_heatflux, rans_const_c, rans_const_sigma,&
     611             rans_mode,                                                        &
    604612             rayleigh_damping_factor, rayleigh_damping_height,                 &
    605613             recycling_width, recycling_yshift,                                &
  • palm/trunk/SOURCE/pmc_interface_mod.f90

    • Property svn:mergeinfo deleted
  • palm/trunk/SOURCE/timestep.f90

    r3049 r3083  
    2525! -----------------
    2626! $Id$
     27! limit dt_3d to be at maximum 2*old_dt; define old_dt at beginning of routine
     28! Add km/kh_max
     29!
     30! 3049 2018-05-29 13:52:36Z Giersch
    2731! Error messages revised
    2832!
     
    161165    INTEGER(iwp) ::  j !<
    162166    INTEGER(iwp) ::  k !<
     167    INTEGER(iwp) ::  km_max_ijk(3) = -1  !< index values (i,j,k) of location where km_max occurs
     168    INTEGER(iwp) ::  kh_max_ijk(3) = -1  !< index values (i,j,k) of location where kh_max occurs
    163169
    164170    LOGICAL ::  stop_dt_local !< local switch for controlling the time stepping
     
    173179    REAL(wp) ::  dt_w              !<
    174180    REAL(wp) ::  dt_w_l            !<
     181    REAL(wp) ::  km_max            !< maximum of Km in entire domain
     182    REAL(wp) ::  kh_max            !< maximum of Kh in entire domain
    175183    REAL(wp) ::  u_gtrans_l        !<
    176184    REAL(wp) ::  u_max_l           !<
     
    190198
    191199
    192 
    193200    CALL cpu_log( log_point(12), 'calculate_timestep', 'start' )
     201!
     202!--    Save former time step as reference
     203       old_dt = dt_3d
    194204
    195205!
     
    318328!--    The time step is the minimum of the 3-4 components and the diffusion time
    319329!--    step minus a reduction (cfl_factor) to be on the safe side.
    320 !--    The time step must not exceed the maximum allowed value.
     330!--    The time step must not exceed the maximum allowed value and must not
     331!--    increase by more than a factor of 10.
    321332       dt_3d = cfl_factor * MIN( dt_diff, dt_u, dt_v, dt_w, dt_precipitation )
    322        dt_3d = MIN( dt_3d, dt_max )
     333       dt_3d = MIN( dt_3d, dt_max, 2.0_wp * old_dt )
    323334
    324335!
     
    334345       IF ( dt_3d < ( 0.00001_wp * dt_max ) )  THEN
    335346          stop_dt = .TRUE.
     347
     348!
     349!--       Determine the maxima of the diffusion coefficients, including their
     350!--       grid index positions.
     351          CALL global_min_max( nzb, nzt+1, nysg, nyng, nxlg, nxrg, km, 'abs',  &
     352                               0.0_wp, km_max, km_max_ijk )
     353          CALL global_min_max( nzb, nzt+1, nysg, nyng, nxlg, nxrg, kh, 'abs',  &
     354                               0.0_wp, kh_max, kh_max_ijk )
    336355
    337356          WRITE( message_string, * ) 'Time step has reached minimum limit.',   &
     
    342361               '&dt_w            = ', dt_w, ' s',                              &
    343362               '&dt_diff         = ', dt_diff, ' s',                           &
    344                '&u_max           = ', u_max, ' m/s   k=', u_max_ijk(1),        &
     363               '&u_max           = ', u_max, ' m/s    k=', u_max_ijk(1),       &
    345364               '  j=', u_max_ijk(2), '  i=', u_max_ijk(3),                     &
    346                '&v_max           = ', v_max, ' m/s   k=', v_max_ijk(1),        &
     365               '&v_max           = ', v_max, ' m/s    k=', v_max_ijk(1),       &
    347366               '  j=', v_max_ijk(2), '  i=', v_max_ijk(3),                     &
    348                '&w_max           = ', w_max, ' m/s   k=', w_max_ijk(1),        &
    349                '  j=', w_max_ijk(2), '  i=', w_max_ijk(3)
     367               '&w_max           = ', w_max, ' m/s    k=', w_max_ijk(1),       &
     368               '  j=', w_max_ijk(2), '  i=', w_max_ijk(3),                     &
     369               '&km_max          = ', km_max, ' m2/s2  k=', km_max_ijk(1),     &
     370               '  j=', km_max_ijk(2), '  i=', km_max_ijk(3),                   &
     371               '&kh_max          = ', kh_max, ' m2/s2  k=', kh_max_ijk(1),     &
     372                '  j=', kh_max_ijk(2), '  i=', kh_max_ijk(3)
    350373          CALL message( 'timestep', 'PA0312', 0, 1, 0, 6, 0 )
    351374!
     
    388411       dt_3d = NINT( dt_3d * 100.0_wp / div ) * div / 100.0_wp
    389412
    390 !
    391 !--    Adjust the time step
    392        old_dt = dt_3d
    393 
    394413    ENDIF
    395414
  • 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 )
  • palm/trunk/SOURCE/vertical_nesting_mod.f90

    r3066 r3083  
    40584058!-- Identical timestep for coarse and fine grids
    40594059          dt_3d = MIN( dtc, dtf )
     4060          !> @fixme setting old_dt might be obsolete at this point
     4061          !>   Due to changes in timestep routine, setting of old_dt might be
     4062          !>   not necessary any more at this point. However, could not be
     4063          !>   tested so far.
     4064          !>   2018-05-18, gronemeier
    40604065          old_dt = dt_3d
    40614066#endif
Note: See TracChangeset for help on using the changeset viewer.