Ignore:
Timestamp:
Oct 22, 2018 7:30:24 PM (6 years ago)
Author:
knoop
Message:

Refactored production_e and production_e_ij (removed code redundancy)

File:
1 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/turbulence_closure_mod.f90

    r3386 r3398  
    2525! -----------------
    2626! $Id$
     27! Refactored production_e and production_e_ij (removed code redundancy)
     28!
     29! 3386 2018-10-19 16:28:22Z gronemeier
    2730! Renamed tcm_prognostic to tcm_prognostic_equations
    2831!
     
    227230    REAL(wp) ::  c_1                !< model constant for RANS mode
    228231    REAL(wp) ::  c_2                !< model constant for RANS mode
    229 !    REAL(wp) ::  c_3                !< model constant for RANS mode
     232    REAL(wp) ::  c_3                !< model constant for RANS mode
    230233    REAL(wp) ::  c_4                !< model constant for RANS mode
    231234    REAL(wp) ::  l_max              !< maximum length scale for Blackadar mixing length
     
    234237
    235238    REAL(wp), DIMENSION(0:4) :: rans_const_c = &       !< model constants for RANS mode (namelist param)
    236        (/ 0.55_wp, 1.44_wp, 1.92_wp, 0.0_wp, 0.0_wp /) !> default values fit for standard-tke-e closure
     239       (/ 0.55_wp, 1.44_wp, 1.92_wp, 1.44_wp, 0.0_wp /) !> default values fit for standard-tke-e closure
    237240
    238241    REAL(wp), DIMENSION(2) :: rans_const_sigma = &     !< model constants for RANS mode, sigma values (sigma_e, sigma_diss) (namelist param)
     
    423426       c_1 = rans_const_c(1)
    424427       c_2 = rans_const_c(2)
    425        !c_3 = rans_const_c(3)   !> @todo clarify how to switch between different models
     428       c_3 = rans_const_c(3)   !> @todo clarify how to switch between different models
    426429       c_4 = rans_const_c(4)
    427430
     
    20942097       IF ( rans_tke_e )  advec = tend
    20952098
    2096        CALL production_e
     2099       CALL production_e( .FALSE. )
    20972100
    20982101!
     
    25472550!> @todo Adjust production term in case of rans_tke_e simulation
    25482551!------------------------------------------------------------------------------!
    2549  SUBROUTINE production_e
     2552 SUBROUTINE production_e( diss_production )
    25502553
    25512554    USE arrays_3d,                                                             &
     
    25682571
    25692572    IMPLICIT NONE
     2573
     2574    LOGICAL :: diss_production
    25702575
    25712576    INTEGER(iwp) ::  i       !< running index x-direction
     
    26002605    REAL(wp), DIMENSION(nzb+1:nzt) ::  dwdy  !< Gradient of w-component in y-direction
    26012606    REAL(wp), DIMENSION(nzb+1:nzt) ::  dwdz  !< Gradient of w-component in z-direction
     2607    REAL(wp), DIMENSION(nzb+1:nzt) ::  tmp_flux  !< temporary flux-array in z-direction
    26022608
    26032609
     
    28242830             flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_0(k,j,i),flag_nr) )
    28252831
    2826              tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def * flag
     2832             IF ( .NOT. diss_production )  THEN
     2833
     2834!--             Compute temdency for TKE-production from shear
     2835                tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def * flag
     2836
     2837             ELSE
     2838
     2839!--             RANS mode: Compute temdency for dissipation-rate-production from shear
     2840                tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def * flag *           &
     2841                              diss(k,j,i)/( e(k,j,i) + 1.0E-20_wp ) * c_1
     2842
     2843             ENDIF
    28272844
    28282845          ENDDO
     
    28442861             DO  i = nxl, nxr
    28452862                DO  j = nys, nyn
     2863
    28462864                   DO  k = nzb+1, nzt
    2847                       tend(k,j,i) = tend(k,j,i) +                              &
    2848                                     kh(k,j,i) * g /                            &
    2849                            MERGE( rho_reference, prho(k,j,i),                  &
    2850                                   use_single_reference_value ) *               &
    2851                                     ( prho(k+1,j,i) - prho(k-1,j,i) ) *        &
    2852                                     dd2zu(k) *                                 &
    2853                                 MERGE( 1.0_wp, 0.0_wp,                         &
    2854                                        BTEST( wall_flags_0(k,j,i), 30 )        &
    2855                                      )                            *            &
    2856                                 MERGE( 1.0_wp, 0.0_wp,                         &
    2857                                        BTEST( wall_flags_0(k,j,i), 9 )         &
    2858                                      )
     2865                      tmp_flux(k) = kh(k,j,i) * ( prho(k+1,j,i) - prho(k-1,j,i) ) * dd2zu(k)
    28592866                   ENDDO
    28602867!
     
    28672874                         DO  m = surf_s, surf_e
    28682875                            k = surf_def_h(l)%k(m)
    2869                             tend(k,j,i) = tend(k,j,i) + g /                    &
    2870                                       MERGE( rho_reference, prho(k,j,i),       &
    2871                                              use_single_reference_value ) *    &
    2872                                       drho_air_zw(k-1) *                       &
    2873                                       surf_def_h(l)%shf(m)
     2876                            tmp_flux(k) = drho_air_zw(k-1) * surf_def_h(l)%shf(m)
    28742877                         ENDDO
    28752878                      ENDDO
    2876 
    28772879                   ENDIF
    28782880
     
    28822884                      DO  m = surf_s, surf_e
    28832885                         k = surf_def_h(2)%k(m)
    2884                          tend(k,j,i) = tend(k,j,i) + g /                       &
    2885                                       MERGE( rho_reference, prho(k,j,i),       &
    2886                                              use_single_reference_value ) *    &
    2887                                       drho_air_zw(k) *                         &
    2888                                       surf_def_h(2)%shf(m)
     2886                         tmp_flux(k) = drho_air_zw(k) * surf_def_h(2)%shf(m)
    28892887                      ENDDO
    28902888                   ENDIF
    28912889
     2890                   IF ( .NOT. diss_production )  THEN
     2891
     2892!--                   Compute temdency for TKE-production from shear
     2893                      DO  k = nzb+1, nzt
     2894                         flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_0(k,j,i),0) )
     2895                         tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g / &
     2896                                       MERGE( rho_reference, prho(k,j,i),       &
     2897                                              use_single_reference_value ) )
     2898                      ENDDO
     2899
     2900                   ELSE
     2901
     2902!--                   RANS mode: Compute temdency for dissipation-rate-production from shear
     2903                      DO  k = nzb+1, nzt
     2904                         flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_0(k,j,i),0) )
     2905                         tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g / &
     2906                                       MERGE( rho_reference, prho(k,j,i),       &
     2907                                              use_single_reference_value ) ) *  &
     2908                                       diss(k,j,i)/( e(k,j,i) + 1.0E-20_wp ) *  &
     2909                                       c_3
     2910                      ENDDO
     2911
     2912                   ENDIF
     2913
    28922914                ENDDO
    28932915             ENDDO
     
    28992921
    29002922                   DO  k = nzb+1, nzt
    2901 !
    2902 !--                   Flag 9 is used to mask top fluxes, flag 30 to mask
    2903 !--                   surface fluxes
    2904                       tend(k,j,i) = tend(k,j,i) -                              &
    2905                                     kh(k,j,i) * g /                            &
    2906                                 MERGE( pt_reference, pt(k,j,i),                &
    2907                                         use_single_reference_value ) *         &
    2908                                     ( pt(k+1,j,i) - pt(k-1,j,i) ) *            &
    2909                                     dd2zu(k)                      *            &
    2910                                 MERGE( 1.0_wp, 0.0_wp,                         &
    2911                                        BTEST( wall_flags_0(k,j,i), 30 )        &
    2912                                      )                            *            &
    2913                                 MERGE( 1.0_wp, 0.0_wp,                         &
    2914                                        BTEST( wall_flags_0(k,j,i), 9 )         &
    2915                                      )
     2923                      tmp_flux(k) = -1.0_wp * kh(k,j,i) * ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k)
    29162924                   ENDDO
    29172925
     
    29242932                         DO  m = surf_s, surf_e
    29252933                            k = surf_def_h(l)%k(m)
    2926                             tend(k,j,i) = tend(k,j,i) + g /                    &
    2927                                  MERGE( pt_reference, pt(k,j,i),               &
    2928                                         use_single_reference_value )           &
    2929                                                    * drho_air_zw(k-1)          &
    2930                                                    * surf_def_h(l)%shf(m)
     2934                            tmp_flux(k) = drho_air_zw(k-1) * surf_def_h(l)%shf(m)
    29312935                         ENDDO
    29322936                      ENDDO
     
    29372941                      DO  m = surf_s, surf_e
    29382942                         k = surf_lsm_h%k(m)
    2939                          tend(k,j,i) = tend(k,j,i) + g /                       &
    2940                                  MERGE( pt_reference, pt(k,j,i),               &
    2941                                         use_single_reference_value )           &
    2942                                                    * drho_air_zw(k-1)          &
    2943                                                    * surf_lsm_h%shf(m)
     2943                         tmp_flux(k) = drho_air_zw(k-1) * surf_lsm_h%shf(m)
    29442944                      ENDDO
    29452945!
     
    29492949                      DO  m = surf_s, surf_e
    29502950                         k = surf_usm_h%k(m)
    2951                          tend(k,j,i) = tend(k,j,i) + g /                       &
    2952                                  MERGE( pt_reference, pt(k,j,i),               &
    2953                                         use_single_reference_value )           &
    2954                                                    * drho_air_zw(k-1)          &
    2955                                                    * surf_usm_h%shf(m)
     2951                         tmp_flux(k) = drho_air_zw(k-1) * surf_usm_h%shf(m)
    29562952                      ENDDO
    29572953                   ENDIF
     
    29622958                      DO  m = surf_s, surf_e
    29632959                         k = surf_def_h(2)%k(m)
    2964                          tend(k,j,i) = tend(k,j,i) + g /                       &
    2965                                  MERGE( pt_reference, pt(k,j,i),               &
    2966                                         use_single_reference_value )           &
    2967                                                    * drho_air_zw(k)            &
    2968                                                    * surf_def_h(2)%shf(m)
     2960                         tmp_flux(k) = drho_air_zw(k) * surf_def_h(2)%shf(m)
    29692961                      ENDDO
    29702962                   ENDIF
    29712963
     2964                   IF ( .NOT. diss_production )  THEN
     2965
     2966!--                   Compute temdency for TKE-production from shear
     2967                      DO  k = nzb+1, nzt
     2968                         flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_0(k,j,i),0) )
     2969                         tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g / &
     2970                                       MERGE( pt_reference, pt(k,j,i),          &
     2971                                              use_single_reference_value ) )
     2972                      ENDDO
     2973
     2974                   ELSE
     2975
     2976!--                   RANS mode: Compute temdency for dissipation-rate-production from shear
     2977                      DO  k = nzb+1, nzt
     2978                         flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_0(k,j,i),0) )
     2979                         tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g / &
     2980                                       MERGE( pt_reference, pt(k,j,i),          &
     2981                                              use_single_reference_value ) ) *  &
     2982                                       diss(k,j,i)/( e(k,j,i) + 1.0E-20_wp ) *  &
     2983                                       c_3
     2984                      ENDDO
     2985
     2986                   ENDIF
     2987
    29722988                ENDDO
    29732989             ENDDO
     
    29812997
    29822998                DO  k = nzb+1, nzt
    2983 !
    2984 !--                Flag 9 is used to mask top fluxes, flag 30 to mask
    2985 !--                surface fluxes
     2999
    29863000                   IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets ) THEN
    29873001                      k1 = 1.0_wp + 0.61_wp * q(k,j,i)
    29883002                      k2 = 0.61_wp * pt(k,j,i)
    2989                       tend(k,j,i) = tend(k,j,i) - kh(k,j,i) *                  &
    2990                                       g /                                      &
    2991                                  MERGE( vpt_reference, vpt(k,j,i),             &
    2992                                         use_single_reference_value ) *         &
    2993                                       ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) +     &
    2994                                         k2 * ( q(k+1,j,i) - q(k-1,j,i) )       &
    2995                                       ) * dd2zu(k) *                           &
    2996                                    MERGE( 1.0_wp, 0.0_wp,                      &
    2997                                           BTEST( wall_flags_0(k,j,i), 30 )     &
    2998                                         )          *                           &
    2999                                    MERGE( 1.0_wp, 0.0_wp,                      &
    3000                                           BTEST( wall_flags_0(k,j,i), 9 )      &
    3001                                         )
     3003                      tmp_flux(k) = -1.0_wp * kh(k,j,i) *                      &
     3004                                      ( k1 * ( pt(k+1,j,i) - pt(k-1,j,i) ) +   &
     3005                                        k2 * ( q(k+1,j,i)  - q(k-1,j,i) )      &
     3006                                      ) * dd2zu(k)
    30023007                   ELSE IF ( bulk_cloud_model )  THEN
    30033008                      IF ( ql(k,j,i) == 0.0_wp )  THEN
     
    30143019                         k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp )
    30153020                      ENDIF
    3016                       tend(k,j,i) = tend(k,j,i) - kh(k,j,i) *                  &
    3017                                       g /                                      &
    3018                                  MERGE( vpt_reference, vpt(k,j,i),             &
    3019                                         use_single_reference_value ) *         &
    3020                                       ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) +     &
    3021                                         k2 * ( q(k+1,j,i) - q(k-1,j,i) )       &
    3022                                       ) * dd2zu(k) *                           &
    3023                                    MERGE( 1.0_wp, 0.0_wp,                      &
    3024                                           BTEST( wall_flags_0(k,j,i), 30 )     &
    3025                                         )          *                           &
    3026                                    MERGE( 1.0_wp, 0.0_wp,                      &
    3027                                           BTEST( wall_flags_0(k,j,i), 9 )      &
    3028                                         )
     3021                      tmp_flux(k) = -1.0_wp * kh(k,j,i) *                      &
     3022                                      ( k1 * ( pt(k+1,j,i) - pt(k-1,j,i) ) +   &
     3023                                        k2 * ( q(k+1,j,i)  - q(k-1,j,i) )      &
     3024                                      ) * dd2zu(k)
    30293025                   ELSE IF ( cloud_droplets )  THEN
    30303026                      k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
    30313027                      k2 = 0.61_wp * pt(k,j,i)
    3032                       tend(k,j,i) = tend(k,j,i) -                              &
    3033                                     kh(k,j,i) * g /                            &
    3034                                  MERGE( vpt_reference, vpt(k,j,i),             &
    3035                                         use_single_reference_value ) *         &
    3036                                     ( k1 * ( pt(k+1,j,i)- pt(k-1,j,i) ) +      &
    3037                                       k2 * ( q(k+1,j,i) -  q(k-1,j,i) ) -      &
    3038                                       pt(k,j,i) * ( ql(k+1,j,i) -              &
    3039                                       ql(k-1,j,i) ) ) * dd2zu(k) *             &
    3040                                    MERGE( 1.0_wp, 0.0_wp,                      &
    3041                                           BTEST( wall_flags_0(k,j,i), 30 )     &
    3042                                         )                        *             &
    3043                                    MERGE( 1.0_wp, 0.0_wp,                      &
    3044                                           BTEST( wall_flags_0(k,j,i), 9 )      &
    3045                                         )
     3028                      tmp_flux(k) = -1.0_wp * kh(k,j,i) * &
     3029                                      ( k1 * ( pt(k+1,j,i) - pt(k-1,j,i) ) +   &
     3030                                        k2 * ( q(k+1,j,i)  - q(k-1,j,i) ) -    &
     3031                                        pt(k,j,i) * ( ql(k+1,j,i) -            &
     3032                                        ql(k-1,j,i) ) ) * dd2zu(k)
    30463033                   ENDIF
    30473034
     
    30793066                         ENDIF
    30803067
    3081                          tend(k,j,i) = tend(k,j,i) + g /                       &
    3082                                  MERGE( vpt_reference, vpt(k,j,i),             &
    3083                                         use_single_reference_value ) *         &
    3084                                             ( k1 * surf_def_h(l)%shf(m) +      &
    3085                                               k2 * surf_def_h(l)%qsws(m)       &
    3086                                             ) * drho_air_zw(k-1)
     3068                         tmp_flux(k) = ( k1 * surf_def_h(l)%shf(m) +           &
     3069                                         k2 * surf_def_h(l)%qsws(m)            &
     3070                                       ) * drho_air_zw(k-1)
    30873071                      ENDDO
    30883072                   ENDDO
     
    31163100                      ENDIF
    31173101
    3118                       tend(k,j,i) = tend(k,j,i) + g /                          &
    3119                                  MERGE( vpt_reference, vpt(k,j,i),             &
    3120                                         use_single_reference_value ) *         &
    3121                                             ( k1 * surf_lsm_h%shf(m) +         &
    3122                                               k2 * surf_lsm_h%qsws(m)          &
    3123                                             ) * drho_air_zw(k-1)
     3102                      tmp_flux(k) = ( k1 * surf_lsm_h%shf(m) +                 &
     3103                                      k2 * surf_lsm_h%qsws(m)                  &
     3104                                    ) * drho_air_zw(k-1)
    31243105                   ENDDO
    31253106!
     
    31523133                      ENDIF
    31533134
    3154                       tend(k,j,i) = tend(k,j,i) + g /                          &
    3155                                  MERGE( vpt_reference, vpt(k,j,i),             &
    3156                                         use_single_reference_value ) *         &
    3157                                             ( k1 * surf_usm_h%shf(m) +         &
    3158                                               k2 * surf_usm_h%qsws(m)          &
    3159                                             ) * drho_air_zw(k-1)
     3135                      tmp_flux(k) = ( k1 * surf_usm_h%shf(m) +                 &
     3136                                      k2 * surf_usm_h%qsws(m)                  &
     3137                                    ) * drho_air_zw(k-1)
    31603138                   ENDDO
    31613139
     
    31913169                      ENDIF
    31923170
    3193                       tend(k,j,i) = tend(k,j,i) + g /                          &
    3194                                  MERGE( vpt_reference, vpt(k,j,i),             &
    3195                                         use_single_reference_value ) *         &
    3196                                             ( k1 * surf_def_h(2)%shf(m) +      &
    3197                                               k2 * surf_def_h(2)%qsws(m)       &
    3198                                             ) * drho_air_zw(k)
     3171                      tmp_flux(k) = ( k1 * surf_def_h(2)%shf(m) +              &
     3172                                      k2 * surf_def_h(2)%qsws(m)               &
     3173                                    ) * drho_air_zw(k)
    31993174
    32003175                   ENDDO
    32013176
    32023177                ENDIF ! from IF ( use_top_fluxes )  THEN
     3178
     3179                IF ( .NOT. diss_production )  THEN
     3180
     3181!--                Compute temdency for TKE-production from shear
     3182                   DO  k = nzb+1, nzt
     3183                      flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_0(k,j,i),0) )
     3184                      tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g / &
     3185                                    MERGE( vpt_reference, vpt(k,j,i),          &
     3186                                           use_single_reference_value ) )
     3187                   ENDDO
     3188
     3189                ELSE
     3190
     3191!--                RANS mode: Compute temdency for dissipation-rate-production from shear
     3192                   DO  k = nzb+1, nzt
     3193                      flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_0(k,j,i),0) )
     3194                      tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g /   &
     3195                                    MERGE( vpt_reference, vpt(k,j,i),          &
     3196                                           use_single_reference_value ) ) *    &
     3197                                    diss(k,j,i)/( e(k,j,i) + 1.0E-20_wp ) *    &
     3198                                    c_3
     3199                   ENDDO
     3200
     3201                ENDIF
    32033202
    32043203             ENDDO
     
    32273226
    32283227    USE control_parameters,                                                    &
    3229         ONLY:  cloud_droplets, constant_flux_layer, neutral,                &
     3228        ONLY:  cloud_droplets, constant_flux_layer, neutral,                   &
    32303229               rho_reference, use_single_reference_value, use_surface_fluxes,  &
    32313230               use_top_fluxes
     
    32393238    USE surface_mod,                                                           &
    32403239        ONLY :  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h,    &
    3241                 surf_usm_v 
     3240                surf_usm_v
    32423241
    32433242    IMPLICIT NONE
     
    32523251    INTEGER(iwp) ::  surf_e  !< end index of surface elements at given i-j position
    32533252    INTEGER(iwp) ::  surf_s  !< start index of surface elements at given i-j position
     3253    INTEGER(iwp) ::  flag_nr !< number of masking flag
    32543254
    32553255    REAL(wp)     ::  def         !<
     
    32603260    REAL(wp)     ::  theta       !<
    32613261    REAL(wp)     ::  temp        !<
    3262     REAL(wp)     ::  sign_dir    !< sign of wall-tke flux, depending on wall orientation 
     3262    REAL(wp)     ::  sign_dir    !< sign of wall-tke flux, depending on wall orientation
    32633263    REAL(wp)     ::  usvs        !< momentum flux u"v"
    32643264    REAL(wp)     ::  vsus        !< momentum flux v"u"
     
    32663266    REAL(wp)     ::  wsvs        !< momentum flux w"v"
    32673267
    3268 
    3269     REAL(wp), DIMENSION(nzb+1:nzt)  ::  dudx        !< Gradient of u-component in x-direction
    3270     REAL(wp), DIMENSION(nzb+1:nzt)  ::  dudy        !< Gradient of u-component in y-direction
    3271     REAL(wp), DIMENSION(nzb+1:nzt)  ::  dudz        !< Gradient of u-component in z-direction
    3272     REAL(wp), DIMENSION(nzb+1:nzt)  ::  dvdx        !< Gradient of v-component in x-direction
    3273     REAL(wp), DIMENSION(nzb+1:nzt)  ::  dvdy        !< Gradient of v-component in y-direction
    3274     REAL(wp), DIMENSION(nzb+1:nzt)  ::  dvdz        !< Gradient of v-component in z-direction
    3275     REAL(wp), DIMENSION(nzb+1:nzt)  ::  dwdx        !< Gradient of w-component in x-direction
    3276     REAL(wp), DIMENSION(nzb+1:nzt)  ::  dwdy        !< Gradient of w-component in y-direction
    3277     REAL(wp), DIMENSION(nzb+1:nzt)  ::  dwdz        !< Gradient of w-component in z-direction
    3278     REAL(wp), DIMENSION(nzb+1:nzt)  ::  tend_temp   !< temporal tendency
     3268    REAL(wp), DIMENSION(nzb+1:nzt) ::  dudx  !< Gradient of u-component in x-direction
     3269    REAL(wp), DIMENSION(nzb+1:nzt) ::  dudy  !< Gradient of u-component in y-direction
     3270    REAL(wp), DIMENSION(nzb+1:nzt) ::  dudz  !< Gradient of u-component in z-direction
     3271    REAL(wp), DIMENSION(nzb+1:nzt) ::  dvdx  !< Gradient of v-component in x-direction
     3272    REAL(wp), DIMENSION(nzb+1:nzt) ::  dvdy  !< Gradient of v-component in y-direction
     3273    REAL(wp), DIMENSION(nzb+1:nzt) ::  dvdz  !< Gradient of v-component in z-direction
     3274    REAL(wp), DIMENSION(nzb+1:nzt) ::  dwdx  !< Gradient of w-component in x-direction
     3275    REAL(wp), DIMENSION(nzb+1:nzt) ::  dwdy  !< Gradient of w-component in y-direction
     3276    REAL(wp), DIMENSION(nzb+1:nzt) ::  dwdz  !< Gradient of w-component in z-direction
     3277    REAL(wp), DIMENSION(nzb+1:nzt) ::  tmp_flux  !< temporary flux-array in z-direction
     3278
     3279
     3280
     3281!
     3282!-- Calculate TKE production by shear. Calculate gradients at all grid
     3283!-- points first, gradients at surface-bounded grid points will be
     3284!-- overwritten further below.
     3285    DO  k = nzb+1, nzt
     3286
     3287       dudx(k) =           ( u(k,j,i+1) - u(k,j,i)     ) * ddx
     3288       dudy(k) = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) -                 &
     3289                             u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
     3290       dudz(k) = 0.5_wp  * ( u(k+1,j,i) + u(k+1,j,i+1) -                 &
     3291                             u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
     3292
     3293       dvdx(k) = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) -                 &
     3294                             v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
     3295       dvdy(k) =           ( v(k,j+1,i) - v(k,j,i)     ) * ddy
     3296       dvdz(k) = 0.5_wp  * ( v(k+1,j,i) + v(k+1,j+1,i) -                 &
     3297                               v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
     3298
     3299       dwdx(k) = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) -                 &
     3300                             w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
     3301       dwdy(k) = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) -                 &
     3302                             w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
     3303       dwdz(k) =           ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
     3304
     3305    ENDDO
     3306
     3307    flag_nr = 29
    32793308
    32803309    IF ( constant_flux_layer )  THEN
    3281 !
    3282 !--    Calculate TKE production by shear. Calculate gradients at all grid
    3283 !--    points first, gradients at surface-bounded grid points will be
    3284 !--    overwritten further below.
    3285        DO  k = nzb+1, nzt
    3286 
    3287           dudx(k)  =           ( u(k,j,i+1) - u(k,j,i)     ) * ddx
    3288           dudy(k)  = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) -                   &
    3289                                  u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
    3290           dudz(k)  = 0.5_wp  * ( u(k+1,j,i) + u(k+1,j,i+1) -                   &
    3291                                  u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
    3292 
    3293           dvdx(k)  = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) -                   &
    3294                                  v(k,j,i-1) - v(k,j+1,i-1) ) * ddx     
    3295           dvdy(k)  =           ( v(k,j+1,i) - v(k,j,i)     ) * ddy
    3296           dvdz(k)  = 0.5_wp  * ( v(k+1,j,i) + v(k+1,j+1,i) -                   &
    3297                                  v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
    3298 
    3299           dwdx(k)  = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) -                   &
    3300                                  w(k,j,i-1) - w(k-1,j,i-1) ) * ddx     
    3301           dwdy(k)  = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) -                   &
    3302                                  w(k,j-1,i) - w(k-1,j-1,i) ) * ddy     
    3303           dwdz(k)  =           ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
    3304 
    3305        ENDDO
     3310
     3311       flag_nr = 0
     3312
     3313!--    Position beneath wall
     3314!--    (2) - Will allways be executed.
     3315!--    'bottom and wall: use u_0,v_0 and wall functions'
    33063316!
    33073317!--    Compute gradients at north- and south-facing surfaces.
    3308 !--    Note, no vertical natural surfaces so far.
     3318!--    First, for default surfaces, then for urban surfaces.
     3319!--    Note, so far no natural vertical surfaces implemented
    33093320       DO  l = 0, 1
    3310 !
    3311 !--       Default surfaces
    33123321          surf_s = surf_def_v(l)%start_index(j,i)
    33133322          surf_e = surf_def_v(l)%end_index(j,i)
     
    33173326             wsvs        = surf_def_v(l)%mom_flux_tke(1,m)
    33183327
    3319              km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp               &
     3328             km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp         &
    33203329                             * 0.5_wp * dy
    33213330!
    33223331!--          -1.0 for right-facing wall, 1.0 for left-facing wall
    3323              sign_dir = MERGE( 1.0_wp, -1.0_wp,                                &
    3324                                BTEST( wall_flags_0(k,j-1,i), 0 ) )
     3332             sign_dir = MERGE( 1.0_wp, -1.0_wp,                          &
     3333                               BTEST( wall_flags_0(k,j-1,i), flag_nr ) )
    33253334             dudy(k) = sign_dir * usvs / ( km_neutral + 1E-10_wp )
    33263335             dwdy(k) = sign_dir * wsvs / ( km_neutral + 1E-10_wp )
    3327           ENDDO   
     3336          ENDDO
    33283337!
    33293338!--       Natural surfaces
     
    33353344             wsvs        = surf_lsm_v(l)%mom_flux_tke(1,m)
    33363345
    3337              km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp               &
     3346             km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp         &
    33383347                             * 0.5_wp * dy
    33393348!
    33403349!--          -1.0 for right-facing wall, 1.0 for left-facing wall
    3341              sign_dir = MERGE( 1.0_wp, -1.0_wp,                                &
    3342                                BTEST( wall_flags_0(k,j-1,i), 0 ) )
     3350             sign_dir = MERGE( 1.0_wp, -1.0_wp,                          &
     3351                               BTEST( wall_flags_0(k,j-1,i), flag_nr ) )
    33433352             dudy(k) = sign_dir * usvs / ( km_neutral + 1E-10_wp )
    33443353             dwdy(k) = sign_dir * wsvs / ( km_neutral + 1E-10_wp )
    3345           ENDDO 
     3354          ENDDO
    33463355!
    33473356!--       Urban surfaces
     
    33533362             wsvs        = surf_usm_v(l)%mom_flux_tke(1,m)
    33543363
    3355              km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp               &
     3364             km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp         &
    33563365                             * 0.5_wp * dy
    33573366!
    33583367!--          -1.0 for right-facing wall, 1.0 for left-facing wall
    3359              sign_dir = MERGE( 1.0_wp, -1.0_wp,                                &
    3360                                BTEST( wall_flags_0(k,j-1,i), 0 ) )
     3368             sign_dir = MERGE( 1.0_wp, -1.0_wp,                          &
     3369                               BTEST( wall_flags_0(k,j-1,i), flag_nr ) )
    33613370             dudy(k) = sign_dir * usvs / ( km_neutral + 1E-10_wp )
    33623371             dwdy(k) = sign_dir * wsvs / ( km_neutral + 1E-10_wp )
    3363           ENDDO 
     3372          ENDDO
    33643373       ENDDO
    33653374!
    3366 !--    Compute gradients at east- and west-facing walls
     3375!--          Compute gradients at east- and west-facing walls
    33673376       DO  l = 2, 3
    3368 !
    3369 !--       Default surfaces
    33703377          surf_s = surf_def_v(l)%start_index(j,i)
    33713378          surf_e = surf_def_v(l)%end_index(j,i)
    33723379          DO  m = surf_s, surf_e
    3373              k           = surf_def_v(l)%k(m)
    3374              vsus        = surf_def_v(l)%mom_flux_tke(0,m)
    3375              wsus        = surf_def_v(l)%mom_flux_tke(1,m)
    3376 
    3377              km_neutral = kappa * ( vsus**2 + wsus**2 )**0.25_wp               &
     3380             k     = surf_def_v(l)%k(m)
     3381             vsus  = surf_def_v(l)%mom_flux_tke(0,m)
     3382             wsus  = surf_def_v(l)%mom_flux_tke(1,m)
     3383
     3384             km_neutral = kappa * ( vsus**2 + wsus**2 )**0.25_wp         &
    33783385                                * 0.5_wp * dx
    33793386!
    33803387!--          -1.0 for right-facing wall, 1.0 for left-facing wall
    3381              sign_dir = MERGE( 1.0_wp, -1.0_wp,                                &
    3382                                BTEST( wall_flags_0(k,j,i-1), 0 ) )
     3388             sign_dir = MERGE( 1.0_wp, -1.0_wp,                          &
     3389                               BTEST( wall_flags_0(k,j,i-1), flag_nr ) )
    33833390             dvdx(k) = sign_dir * vsus / ( km_neutral + 1E-10_wp )
    33843391             dwdx(k) = sign_dir * wsus / ( km_neutral + 1E-10_wp )
    3385           ENDDO 
    3386 !
    3387 !--       Natural surfaces
     3392          ENDDO
     3393!
     3394!--             Natural surfaces
    33883395          surf_s = surf_lsm_v(l)%start_index(j,i)
    33893396          surf_e = surf_lsm_v(l)%end_index(j,i)
    33903397          DO  m = surf_s, surf_e
    3391              k           = surf_lsm_v(l)%k(m)
    3392              vsus        = surf_lsm_v(l)%mom_flux_tke(0,m)
    3393              wsus        = surf_lsm_v(l)%mom_flux_tke(1,m)
    3394 
    3395              km_neutral = kappa * ( vsus**2 + wsus**2 )**0.25_wp               &
     3398             k     = surf_lsm_v(l)%k(m)
     3399             vsus  = surf_lsm_v(l)%mom_flux_tke(0,m)
     3400             wsus  = surf_lsm_v(l)%mom_flux_tke(1,m)
     3401
     3402             km_neutral = kappa * ( vsus**2 + wsus**2 )**0.25_wp         &
    33963403                                * 0.5_wp * dx
    33973404!
    33983405!--          -1.0 for right-facing wall, 1.0 for left-facing wall
    3399              sign_dir = MERGE( 1.0_wp, -1.0_wp,                                &
    3400                                BTEST( wall_flags_0(k,j,i-1), 0 ) )
     3406             sign_dir = MERGE( 1.0_wp, -1.0_wp,                          &
     3407                               BTEST( wall_flags_0(k,j,i-1), flag_nr ) )
    34013408             dvdx(k) = sign_dir * vsus / ( km_neutral + 1E-10_wp )
    34023409             dwdx(k) = sign_dir * wsus / ( km_neutral + 1E-10_wp )
    3403           ENDDO 
    3404 !
    3405 !--       Urban surfaces
     3410          ENDDO
     3411!
     3412!--             Urban surfaces
    34063413          surf_s = surf_usm_v(l)%start_index(j,i)
    34073414          surf_e = surf_usm_v(l)%end_index(j,i)
    34083415          DO  m = surf_s, surf_e
    3409              k           = surf_usm_v(l)%k(m)
    3410              vsus        = surf_usm_v(l)%mom_flux_tke(0,m)
    3411              wsus        = surf_usm_v(l)%mom_flux_tke(1,m)
    3412 
    3413              km_neutral = kappa * ( vsus**2 + wsus**2 )**0.25_wp               &
     3416             k     = surf_usm_v(l)%k(m)
     3417             vsus  = surf_usm_v(l)%mom_flux_tke(0,m)
     3418             wsus  = surf_usm_v(l)%mom_flux_tke(1,m)
     3419
     3420             km_neutral = kappa * ( vsus**2 + wsus**2 )**0.25_wp         &
    34143421                                * 0.5_wp * dx
    34153422!
    34163423!--          -1.0 for right-facing wall, 1.0 for left-facing wall
    3417              sign_dir = MERGE( 1.0_wp, -1.0_wp,                                &
    3418                                BTEST( wall_flags_0(k,j,i-1), 0 ) )
     3424             sign_dir = MERGE( 1.0_wp, -1.0_wp,                          &
     3425                               BTEST( wall_flags_0(k,j,i-1), flag_nr ) )
    34193426             dvdx(k) = sign_dir * vsus / ( km_neutral + 1E-10_wp )
    34203427             dwdx(k) = sign_dir * wsus / ( km_neutral + 1E-10_wp )
    3421           ENDDO 
     3428          ENDDO
    34223429       ENDDO
    34233430!
    3424 !--    Compute gradients at upward-facing walls, first for
    3425 !--    non-natural default surfaces
     3431!--          Compute gradients at upward-facing surfaces
    34263432       surf_s = surf_def_h(0)%start_index(j,i)
    34273433       surf_e = surf_def_h(0)%end_index(j,i)
     
    34293435          k = surf_def_h(0)%k(m)
    34303436!
    3431 !--       Please note, actually, an interpolation of u_0 and v_0 
    3432 !--       onto the grid center would be required. However, this 
     3437!--       Please note, actually, an interpolation of u_0 and v_0
     3438!--       onto the grid center would be required. However, this
    34333439!--       would require several data transfers between 2D-grid and
    3434 !--       wall type. The effect of this missing interpolation is 
     3440!--       wall type. The effect of this missing interpolation is
    34353441!--       negligible. (See also production_e_init).
    3436           dudz(k)     = ( u(k+1,j,i) - surf_def_h(0)%u_0(m) ) * dd2zu(k)   
    3437           dvdz(k)     = ( v(k+1,j,i) - surf_def_h(0)%v_0(m) ) * dd2zu(k)
     3442          dudz(k) = ( u(k+1,j,i) - surf_def_h(0)%u_0(m) ) * dd2zu(k)
     3443          dvdz(k) = ( v(k+1,j,i) - surf_def_h(0)%v_0(m) ) * dd2zu(k)
    34383444
    34393445       ENDDO
     
    34453451          k = surf_lsm_h%k(m)
    34463452
    3447           dudz(k)     = ( u(k+1,j,i) - surf_lsm_h%u_0(m) ) * dd2zu(k)   
    3448           dvdz(k)     = ( v(k+1,j,i) - surf_lsm_h%v_0(m) ) * dd2zu(k)
     3453          dudz(k) = ( u(k+1,j,i) - surf_lsm_h%u_0(m) ) * dd2zu(k)
     3454          dvdz(k) = ( v(k+1,j,i) - surf_lsm_h%v_0(m) ) * dd2zu(k)
     3455
    34493456       ENDDO
    34503457!
     
    34553462          k = surf_usm_h%k(m)
    34563463
    3457           dudz(k)     = ( u(k+1,j,i) - surf_usm_h%u_0(m) ) * dd2zu(k)   
    3458           dvdz(k)     = ( v(k+1,j,i) - surf_usm_h%v_0(m) ) * dd2zu(k)
     3464          dudz(k) = ( u(k+1,j,i) - surf_usm_h%u_0(m) ) * dd2zu(k)
     3465          dvdz(k) = ( v(k+1,j,i) - surf_usm_h%v_0(m) ) * dd2zu(k)
     3466
    34593467       ENDDO
    34603468!
    3461 !--    Compute gradients at downward-facing walls, only for 
     3469!--    Compute gradients at downward-facing walls, only for
    34623470!--    non-natural default surfaces
    34633471       surf_s = surf_def_h(1)%start_index(j,i)
     
    34663474          k = surf_def_h(1)%k(m)
    34673475
    3468           dudz(k)     = ( surf_def_h(1)%u_0(m) - u(k-1,j,i) ) * dd2zu(k)   
    3469           dvdz(k)     = ( surf_def_h(1)%v_0(m) - v(k-1,j,i) ) * dd2zu(k)
     3476          dudz(k) = ( surf_def_h(1)%u_0(m) - u(k-1,j,i) ) * dd2zu(k)
     3477          dvdz(k) = ( surf_def_h(1)%v_0(m) - v(k-1,j,i) ) * dd2zu(k)
    34703478
    34713479       ENDDO
    34723480
    3473 !        IF ( .NOT. rans_tke_e )  THEN
    3474 
    3475        DO  k = nzb+1, nzt
    3476 
    3477           def = 2.0_wp * ( dudx(k)**2 + dvdy(k)**2 + dwdz(k)**2 ) +         &
    3478                            dudy(k)**2 + dvdx(k)**2 + dwdx(k)**2   +         &
    3479                            dwdy(k)**2 + dudz(k)**2 + dvdz(k)**2   +         &
    3480                 2.0_wp * ( dvdx(k)*dudy(k) +                                &
    3481                            dwdx(k)*dudz(k) +                                &
    3482                            dwdy(k)*dvdz(k) )
    3483 
    3484           IF ( def < 0.0_wp )  def = 0.0_wp
    3485 
    3486           flag  = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
    3487 
    3488           tend_temp(k) = km(k,j,i) * def * flag
    3489 
    3490        ENDDO
    3491 
    3492 !        ELSE
    3493 !
    3494 !           DO  k = nzb+1, nzt
    3495 ! !
    3496 ! !--          Production term according to Kato and Launder (1993)
    3497 !              def = SQRT( ( dudy(k)**2 + dvdz(k)**2 + dwdx(k)**2 +              &
    3498 !                            dudz(k)**2 + dvdx(k)**2 + dwdy(k)**2 +              &
    3499 !                            2.0_wp * ( dudy(k) * dvdx(k) +                      &
    3500 !                                       dvdz(k) * dwdy(k) +                      &
    3501 !                                       dwdx(k) * dudz(k) )       )              &
    3502 !                        * ( dudy(k)**2 + dvdz(k)**2 + dwdx(k)**2 +              &
    3503 !                            dudz(k)**2 + dvdx(k)**2 + dwdy(k)**2 -              &
    3504 !                            2.0_wp * ( dudy(k) * dvdx(k) +                      &
    3505 !                                       dvdz(k) * dwdy(k) +                      &
    3506 !                                       dwdx(k) * dudz(k) )       )              &
    3507 !                        )
    3508 !
    3509 !              IF ( def < 0.0_wp )  def = 0.0_wp
    3510 !
    3511 !              flag  = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
    3512 !
    3513 !              tend_temp(k) = km(k,j,i) * def * flag
    3514 !
    3515 !           ENDDO
    3516 !
    3517 !        ENDIF
    3518 
    3519     ELSE  ! not constant_flux_layer
    3520 
    3521 !        IF ( .NOT. rans_tke_e )  THEN
    3522 !
    3523 !--       Calculate TKE production by shear. Here, no additional
    3524 !--       wall-bounded code is considered.
    3525 !--       Why?
    3526           DO  k = nzb+1, nzt
    3527 
    3528              dudx(k)  =           ( u(k,j,i+1) - u(k,j,i)     ) * ddx
    3529              dudy(k)  = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) -                &
    3530                                     u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
    3531              dudz(k)  = 0.5_wp  * ( u(k+1,j,i) + u(k+1,j,i+1) -                &
    3532                                     u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
    3533 
    3534              dvdx(k)  = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) -                &
    3535                                     v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
    3536              dvdy(k)  =           ( v(k,j+1,i) - v(k,j,i)     ) * ddy
    3537              dvdz(k)  = 0.5_wp  * ( v(k+1,j,i) + v(k+1,j+1,i) -                &
    3538                                     v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
    3539 
    3540              dwdx(k)  = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) -                &
    3541                                     w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
    3542              dwdy(k)  = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) -                &
    3543                                     w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
    3544              dwdz(k)  =           ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
    3545 
    3546              def = 2.0_wp * ( dudx(k)**2 + dvdy(k)**2 + dwdz(k)**2 ) +         &
    3547                               dudy(k)**2 + dvdx(k)**2 + dwdx(k)**2   +         &
    3548                               dwdy(k)**2 + dudz(k)**2 + dvdz(k)**2   +         &
    3549                    2.0_wp * ( dvdx(k)*dudy(k) +                                &
    3550                               dwdx(k)*dudz(k) +                                &
    3551                               dwdy(k)*dvdz(k) )
    3552 
    3553              IF ( def < 0.0_wp )  def = 0.0_wp
    3554 
    3555              flag        = MERGE( 1.0_wp, 0.0_wp,                              &
    3556                                   BTEST( wall_flags_0(k,j,i), 29 ) )
    3557              tend_temp(k) = km(k,j,i) * def * flag
    3558 
    3559           ENDDO
    3560 
    3561 !        ELSE
    3562 !
    3563 !           DO  k = nzb+1, nzt
    3564 !
    3565 !              dudx(k)  =           ( u(k,j,i+1) - u(k,j,i)     ) * ddx
    3566 !              dudy(k)  = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) -                &
    3567 !                                     u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
    3568 !              dudz(k)  = 0.5_wp  * ( u(k+1,j,i) + u(k+1,j,i+1) -                &
    3569 !                                     u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
    3570 !
    3571 !              dvdx(k)  = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) -                &
    3572 !                                     v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
    3573 !              dvdy(k)  =           ( v(k,j+1,i) - v(k,j,i)     ) * ddy
    3574 !              dvdz(k)  = 0.5_wp  * ( v(k+1,j,i) + v(k+1,j+1,i) -                &
    3575 !                                     v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
    3576 !
    3577 !              dwdx(k)  = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) -                &
    3578 !                                     w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
    3579 !              dwdy(k)  = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) -                &
    3580 !                                     w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
    3581 !              dwdz(k)  =           ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
    3582 ! !
    3583 ! !--          Production term according to Kato and Launder (1993)
    3584 !              def = SQRT( ( dudy(k)**2 + dvdz(k)**2 + dwdx(k)**2 +              &
    3585 !                            dudz(k)**2 + dvdx(k)**2 + dwdy(k)**2 +              &
    3586 !                            2.0_wp * ( dudy(k) * dvdx(k) +                      &
    3587 !                                       dvdz(k) * dwdy(k) +                      &
    3588 !                                       dwdx(k) * dudz(k) )       )              &
    3589 !                        * ( dudy(k)**2 + dvdz(k)**2 + dwdx(k)**2 +              &
    3590 !                            dudz(k)**2 + dvdx(k)**2 + dwdy(k)**2 -              &
    3591 !                            2.0_wp * ( dudy(k) * dvdx(k) +                      &
    3592 !                                       dvdz(k) * dwdy(k) +                      &
    3593 !                                       dwdx(k) * dudz(k) )       )              &
    3594 !                        )
    3595 !
    3596 !              IF ( def < 0.0_wp )  def = 0.0_wp
    3597 !
    3598 !              flag        = MERGE( 1.0_wp, 0.0_wp,                              &
    3599 !                                   BTEST( wall_flags_0(k,j,i), 29 ) )
    3600 !              tend_temp(k) = km(k,j,i) * def * flag
    3601 !
    3602 !           ENDDO
    3603 !
    3604 !        ENDIF
    3605 
    36063481    ENDIF
    36073482
    3608     IF ( .NOT. diss_production )  THEN
    3609 !
    3610 !--    Production term in case of TKE production
    3611        DO  k = nzb+1, nzt
    3612           tend(k,j,i) = tend(k,j,i) + tend_temp(k)
    3613        ENDDO
    3614     ELSE
    3615 !
    3616 !--    Production term in case of dissipation-rate production (rans_tke_e)
    3617        DO  k = nzb+1, nzt
    3618 
    3619           ! Standard TKE-e closure
    3620           tend(k,j,i) = tend(k,j,i) + tend_temp(k) * diss(k,j,i)               &
    3621                                     /( e(k,j,i) + 1.0E-20_wp )                 &
    3622                                     * c_1
    3623 !           ! Production according to Koblitz (2013)
    3624 !           tend(k,j,i) = tend(k,j,i) + tend_temp(k) * diss(k,j,i)               &
    3625 !                                     /( e(k,j,i) + 1.0E-20_wp )                 &
    3626 !                                     * ( c_1 + ( c_2 - c_1 )                    &
    3627 !                                             * l_wall(k,j,i) / l_max )
    3628 !           ! Production according to Detering and Etling (1985)
    3629 !           !> @todo us is not correct if there are vertical walls
    3630 !           tend(k,j,i) = tend(k,j,i) + tend_temp(k) * SQRT(e(k,j,i))            &
    3631 !                                     * c_1 * c_0**3 / c_4 * f                   &
    3632 !              / surf_def_h(0)%us(surf_def_h(0)%start_index(j,i))
    3633        ENDDO
    3634     ENDIF
     3483    DO  k = nzb+1, nzt
     3484
     3485       def = 2.0_wp * ( dudx(k)**2 + dvdy(k)**2 + dwdz(k)**2 ) +         &
     3486                        dudy(k)**2 + dvdx(k)**2 + dwdx(k)**2 +           &
     3487                        dwdy(k)**2 + dudz(k)**2 + dvdz(k)**2 +           &
     3488             2.0_wp * ( dvdx(k)*dudy(k) + dwdx(k)*dudz(k) +              &
     3489                        dwdy(k)*dvdz(k) )
     3490
     3491       IF ( def < 0.0_wp )  def = 0.0_wp
     3492
     3493       flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_0(k,j,i),flag_nr) )
     3494
     3495       IF ( .NOT. diss_production )  THEN
     3496
     3497!--       Compute temdency for TKE-production from shear
     3498          tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def * flag
     3499
     3500       ELSE
     3501
     3502!--       RANS mode: Compute temdency for dissipation-rate-production from shear
     3503          tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def * flag *           &
     3504                        diss(k,j,i)/( e(k,j,i) + 1.0E-20_wp ) * c_1
     3505
     3506       ENDIF
     3507
     3508    ENDDO
     3509
    36353510
    36363511!
     
    36423517          IF ( ocean_mode )  THEN
    36433518!
    3644 !--          So far in the ocean no special treatment of density flux in
    3645 !--          the bottom and top surface layer
     3519!--          So far in the ocean no special treatment of density flux
     3520!--          in the bottom and top surface layer
     3521
    36463522             DO  k = nzb+1, nzt
    3647  
    3648                 tend(k,j,i) = tend(k,j,i) +                                    &
    3649                               kh(k,j,i) * g /                                  &
    3650                               MERGE( rho_reference, prho(k,j,i),               &
    3651                                      use_single_reference_value ) *            &
    3652                               ( prho(k+1,j,i) - prho(k-1,j,i) ) *              &
    3653                               dd2zu(k) *                                       &
    3654                                 MERGE( 1.0_wp, 0.0_wp,                         &
    3655                                        BTEST( wall_flags_0(k,j,i), 30 )        &
    3656                                      ) *                                       &
    3657                                 MERGE( 1.0_wp, 0.0_wp,                         &
    3658                                        BTEST( wall_flags_0(k,j,i), 9 )         &
    3659                                      )           
     3523                tmp_flux(k) = kh(k,j,i) * ( prho(k+1,j,i) - prho(k-1,j,i) ) * dd2zu(k)
     3524             ENDDO
     3525!
     3526!--          Treatment of near-surface grid points, at up- and down-
     3527!--          ward facing surfaces
     3528             IF ( use_surface_fluxes )  THEN
     3529                DO  l = 0, 1
     3530                   surf_s = surf_def_h(l)%start_index(j,i)
     3531                   surf_e = surf_def_h(l)%end_index(j,i)
     3532                   DO  m = surf_s, surf_e
     3533                      k = surf_def_h(l)%k(m)
     3534                      tmp_flux(k) = drho_air_zw(k-1) * surf_def_h(l)%shf(m)
     3535                   ENDDO
     3536                ENDDO
     3537             ENDIF
     3538
     3539             IF ( use_top_fluxes )  THEN
     3540                surf_s = surf_def_h(2)%start_index(j,i)
     3541                surf_e = surf_def_h(2)%end_index(j,i)
     3542                DO  m = surf_s, surf_e
     3543                   k = surf_def_h(2)%k(m)
     3544                   tmp_flux(k) = drho_air_zw(k) * surf_def_h(2)%shf(m)
     3545                ENDDO
     3546             ENDIF
     3547
     3548             IF ( .NOT. diss_production )  THEN
     3549
     3550!--             Compute temdency for TKE-production from shear
     3551                DO  k = nzb+1, nzt
     3552                   flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_0(k,j,i),0) )
     3553                   tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g / &
     3554                                 MERGE( rho_reference, prho(k,j,i),       &
     3555                                        use_single_reference_value ) )
     3556                ENDDO
     3557
     3558             ELSE
     3559
     3560!--             RANS mode: Compute temdency for dissipation-rate-production from shear
     3561                DO  k = nzb+1, nzt
     3562                   flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_0(k,j,i),0) )
     3563                   tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g / &
     3564                                 MERGE( rho_reference, prho(k,j,i),       &
     3565                                        use_single_reference_value ) ) *  &
     3566                                 diss(k,j,i)/( e(k,j,i) + 1.0E-20_wp ) *  &
     3567                                 c_3
     3568                ENDDO
     3569
     3570             ENDIF
     3571
     3572
     3573          ELSE ! or IF ( .NOT. ocean_mode )  THEN
     3574
     3575             DO  k = nzb+1, nzt
     3576                tmp_flux(k) = -1.0_wp * kh(k,j,i) * ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k)
    36603577             ENDDO
    36613578
     
    36683585                   DO  m = surf_s, surf_e
    36693586                      k = surf_def_h(l)%k(m)
    3670                       tend(k,j,i) = tend(k,j,i) + g /                          &
    3671                                 MERGE( rho_reference, prho(k,j,i),             &
    3672                                        use_single_reference_value ) *          &
    3673                                 drho_air_zw(k-1) *                             &
    3674                                 surf_def_h(l)%shf(m)
     3587                      tmp_flux(k) = drho_air_zw(k-1) * surf_def_h(l)%shf(m)
    36753588                   ENDDO
    36763589                ENDDO
    3677 
     3590!
     3591!--             Natural surfaces
     3592                surf_s = surf_lsm_h%start_index(j,i)
     3593                surf_e = surf_lsm_h%end_index(j,i)
     3594                DO  m = surf_s, surf_e
     3595                   k = surf_lsm_h%k(m)
     3596                   tmp_flux(k) = drho_air_zw(k-1) * surf_lsm_h%shf(m)
     3597                ENDDO
     3598!
     3599!--             Urban surfaces
     3600                surf_s = surf_usm_h%start_index(j,i)
     3601                surf_e = surf_usm_h%end_index(j,i)
     3602                DO  m = surf_s, surf_e
     3603                   k = surf_usm_h%k(m)
     3604                   tmp_flux(k) = drho_air_zw(k-1) * surf_usm_h%shf(m)
     3605                ENDDO
    36783606             ENDIF
    36793607
     
    36833611                DO  m = surf_s, surf_e
    36843612                   k = surf_def_h(2)%k(m)
    3685                    tend(k,j,i) = tend(k,j,i) + g /                             &
    3686                                 MERGE( rho_reference, prho(k,j,i),             &
    3687                                        use_single_reference_value ) *          &
    3688                                 drho_air_zw(k) *                               &
    3689                                 surf_def_h(2)%shf(m)
     3613                   tmp_flux(k) = drho_air_zw(k) * surf_def_h(2)%shf(m)
    36903614                ENDDO
    36913615             ENDIF
    36923616
    3693           ELSE
    3694 
    3695              DO  k = nzb+1, nzt
    3696 !
    3697 !--             Flag 9 is used to mask top fluxes, flag 30 to mask
    3698 !--             surface fluxes
    3699                 tend(k,j,i) = tend(k,j,i) -                                    &
    3700                               kh(k,j,i) * g /                                  &
    3701                                 MERGE( pt_reference, pt(k,j,i),                &
    3702                                        use_single_reference_value ) *          &
    3703                               ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k) *       &
    3704                                 MERGE( 1.0_wp, 0.0_wp,                         &
    3705                                        BTEST( wall_flags_0(k,j,i), 30 )        &
    3706                                      ) *                                       &
    3707                                 MERGE( 1.0_wp, 0.0_wp,                         &
    3708                                        BTEST( wall_flags_0(k,j,i), 9 )         &
    3709                                      )
    3710 
    3711              ENDDO
    3712 
    3713              IF ( use_surface_fluxes )  THEN
    3714 !
    3715 !--             Default surfaces, up- and downward-facing
    3716                 DO  l = 0, 1
    3717                    surf_s = surf_def_h(l)%start_index(j,i)
    3718                    surf_e = surf_def_h(l)%end_index(j,i)
    3719                    DO  m = surf_s, surf_e
    3720                       k = surf_def_h(l)%k(m)
    3721                       tend(k,j,i) = tend(k,j,i) + g /                          &
    3722                                 MERGE( pt_reference, pt(k,j,i),                &
    3723                                        use_single_reference_value ) *          &
    3724                                 drho_air_zw(k-1) *                             &
    3725                                 surf_def_h(l)%shf(m)
    3726                    ENDDO
     3617             IF ( .NOT. diss_production )  THEN
     3618
     3619!--             Compute temdency for TKE-production from shear
     3620                DO  k = nzb+1, nzt
     3621                   flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_0(k,j,i),0) )
     3622                   tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g / &
     3623                                 MERGE( pt_reference, pt(k,j,i),          &
     3624                                        use_single_reference_value ) )
    37273625                ENDDO
    3728 !
    3729 !--             Natural surfaces
    3730                 surf_s = surf_lsm_h%start_index(j,i)
    3731                 surf_e = surf_lsm_h%end_index(j,i)
    3732                 DO  m = surf_s, surf_e
    3733                    k = surf_lsm_h%k(m)
    3734                    tend(k,j,i) = tend(k,j,i) + g /                            &
    3735                                 MERGE( pt_reference, pt(k,j,i),                &
    3736                                        use_single_reference_value ) *          &
    3737                                 drho_air_zw(k-1) *                             &
    3738                                 surf_lsm_h%shf(m)
     3626
     3627             ELSE
     3628
     3629!--             RANS mode: Compute temdency for dissipation-rate-production from shear
     3630                DO  k = nzb+1, nzt
     3631                   flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_0(k,j,i),0) )
     3632                   tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g / &
     3633                                 MERGE( pt_reference, pt(k,j,i),          &
     3634                                        use_single_reference_value ) ) *  &
     3635                                 diss(k,j,i)/( e(k,j,i) + 1.0E-20_wp ) *  &
     3636                                 c_3
    37393637                ENDDO
    3740 !
    3741 !--             Urban surfaces
    3742                 surf_s = surf_usm_h%start_index(j,i)
    3743                 surf_e = surf_usm_h%end_index(j,i)
    3744                 DO  m = surf_s, surf_e
    3745                    k = surf_usm_h%k(m)
    3746                    tend(k,j,i) = tend(k,j,i) + g /                             &
    3747                                 MERGE( pt_reference, pt(k,j,i),                &
    3748                                        use_single_reference_value ) *          &
    3749                                 drho_air_zw(k-1) *                             &
    3750                                 surf_usm_h%shf(m)
    3751                 ENDDO
     3638
    37523639             ENDIF
    37533640
    3754              IF ( use_top_fluxes )  THEN
    3755                 surf_s = surf_def_h(2)%start_index(j,i)
    3756                 surf_e = surf_def_h(2)%end_index(j,i)
    3757                 DO  m = surf_s, surf_e
    3758                    k = surf_def_h(2)%k(m)
    3759                    tend(k,j,i) = tend(k,j,i) + g /                             &
    3760                                 MERGE( pt_reference, pt(k,j,i),                &
    3761                                        use_single_reference_value ) *          &
    3762                                 drho_air_zw(k) *                               &
    3763                                 surf_def_h(2)%shf(m)
    3764                 ENDDO
    3765              ENDIF
    3766 
    3767           ENDIF
    3768 
    3769        ELSE
     3641          ENDIF ! from IF ( .NOT. ocean_mode )
     3642
     3643       ELSE ! or IF ( humidity )  THEN
    37703644
    37713645          DO  k = nzb+1, nzt
    3772 !
    3773 !--          Flag 9 is used to mask top fluxes, flag 30 to mask surface fluxes
    3774              IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets )  THEN
     3646
     3647             IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets ) THEN
    37753648                k1 = 1.0_wp + 0.61_wp * q(k,j,i)
    37763649                k2 = 0.61_wp * pt(k,j,i)
    3777                 tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g /                    &
    3778                                 MERGE( vpt_reference, vpt(k,j,i),              &
    3779                                        use_single_reference_value ) *          &
    3780                                       ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) +     &
    3781                                         k2 * ( q(k+1,j,i) - q(k-1,j,i) )       &
    3782                                       ) * dd2zu(k) *                           &
    3783                                    MERGE( 1.0_wp, 0.0_wp,                      &
    3784                                           BTEST( wall_flags_0(k,j,i), 30 )     &
    3785                                         )          *                           &
    3786                                    MERGE( 1.0_wp, 0.0_wp,                      &
    3787                                           BTEST( wall_flags_0(k,j,i), 9 )      &
    3788                                         )
    3789 
     3650                tmp_flux(k) = -1.0_wp * kh(k,j,i) *                      &
     3651                                ( k1 * ( pt(k+1,j,i) - pt(k-1,j,i) ) +   &
     3652                                  k2 * ( q(k+1,j,i)  - q(k-1,j,i) )      &
     3653                                ) * dd2zu(k)
    37903654             ELSE IF ( bulk_cloud_model )  THEN
    37913655                IF ( ql(k,j,i) == 0.0_wp )  THEN
     
    37953659                   theta = pt(k,j,i) + d_exner(k) * lv_d_cp * ql(k,j,i)
    37963660                   temp  = theta * exner(k)
    3797                    k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *                        &
    3798                                  ( q(k,j,i) - ql(k,j,i) ) *                    &
    3799                         ( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) /              &
    3800                         ( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp *               &
     3661                   k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *                  &
     3662                                 ( q(k,j,i) - ql(k,j,i) ) *              &
     3663                        ( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) /        &
     3664                        ( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp *         &
    38013665                        ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
    38023666                   k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp )
    38033667                ENDIF
    3804                 tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g /                    &
    3805                                 MERGE( vpt_reference, vpt(k,j,i),              &
    3806                                        use_single_reference_value ) *          &
    3807                                       ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) +     &
    3808                                         k2 * ( q(k+1,j,i) - q(k-1,j,i) )       &
    3809                                       ) * dd2zu(k) *                           &
    3810                                    MERGE( 1.0_wp, 0.0_wp,                      &
    3811                                           BTEST( wall_flags_0(k,j,i), 30 )     &
    3812                                         )          *                           &
    3813                                    MERGE( 1.0_wp, 0.0_wp,                      &
    3814                                           BTEST( wall_flags_0(k,j,i), 9 )      &
    3815                                         )
     3668                tmp_flux(k) = -1.0_wp * kh(k,j,i) *                      &
     3669                                ( k1 * ( pt(k+1,j,i) - pt(k-1,j,i) ) +   &
     3670                                  k2 * ( q(k+1,j,i)  - q(k-1,j,i) )      &
     3671                                ) * dd2zu(k)
    38163672             ELSE IF ( cloud_droplets )  THEN
    38173673                k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
    38183674                k2 = 0.61_wp * pt(k,j,i)
    3819                 tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g /                    &
    3820                                 MERGE( vpt_reference, vpt(k,j,i),              &
    3821                                        use_single_reference_value ) *          &
    3822                                   ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) +         &
    3823                                     k2 * ( q(k+1,j,i) - q(k-1,j,i) ) -         &
    3824                                     pt(k,j,i) * ( ql(k+1,j,i) -                &
    3825                                                   ql(k-1,j,i) ) ) * dd2zu(k)   &
    3826                                  * MERGE( 1.0_wp, 0.0_wp,                      &
    3827                                           BTEST( wall_flags_0(k,j,i), 30 )     &
    3828                                         )                                      &
    3829                                  * MERGE( 1.0_wp, 0.0_wp,                      &
    3830                                           BTEST( wall_flags_0(k,j,i), 9 )      &
    3831                                         )   
     3675                tmp_flux(k) = -1.0_wp * kh(k,j,i) * &
     3676                                ( k1 * ( pt(k+1,j,i) - pt(k-1,j,i) ) +   &
     3677                                  k2 * ( q(k+1,j,i)  - q(k-1,j,i) ) -    &
     3678                                  pt(k,j,i) * ( ql(k+1,j,i) -            &
     3679                                  ql(k-1,j,i) ) ) * dd2zu(k)
    38323680             ENDIF
     3681
    38333682          ENDDO
    38343683
    38353684          IF ( use_surface_fluxes )  THEN
    38363685!
    3837 !--          Treat horizontal default surfaces, up- and downward-facing
     3686!--          Treat horizontal default surfaces
    38383687             DO  l = 0, 1
    38393688                surf_s = surf_def_h(l)%start_index(j,i)
     
    38503699                         k2 = 0.61_wp * pt(k,j,i)
    38513700                      ELSE
    3852                         theta = pt(k,j,i) + d_exner(k) * lv_d_cp * ql(k,j,i)
    3853                         temp  = theta * exner(k)
    3854                         k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *                   &
    3855                                    ( q(k,j,i) - ql(k,j,i) ) *                  &
    3856                           ( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) /            &
    3857                           ( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp *             &
    3858                           ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
    3859                         k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp )
     3701                         theta = pt(k,j,i) + d_exner(k) * lv_d_cp * ql(k,j,i)
     3702                         temp  = theta * exner(k)
     3703                         k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *            &
     3704                                    ( q(k,j,i) - ql(k,j,i) ) *           &
     3705                           ( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) /     &
     3706                           ( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp *      &
     3707                           ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
     3708                         k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp )
    38603709                      ENDIF
    38613710                   ELSE IF ( cloud_droplets )  THEN
     
    38643713                   ENDIF
    38653714
    3866                    tend(k,j,i) = tend(k,j,i) + g /                             &
    3867                                 MERGE( vpt_reference, vpt(k,j,i),              &
    3868                                        use_single_reference_value ) *          &
    3869                                       ( k1 * surf_def_h(l)%shf(m) +            &
    3870                                         k2 * surf_def_h(l)%qsws(m)             &
    3871                                       ) * drho_air_zw(k-1)
     3715                   tmp_flux(k) = ( k1 * surf_def_h(l)%shf(m) +           &
     3716                                   k2 * surf_def_h(l)%qsws(m)            &
     3717                                 ) * drho_air_zw(k-1)
    38723718                ENDDO
    38733719             ENDDO
     
    38803726
    38813727                IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets ) THEN
    3882                     k1 = 1.0_wp + 0.61_wp * q(k,j,i)
    3883                     k2 = 0.61_wp * pt(k,j,i)
    3884                 ELSE IF ( bulk_cloud_model )  THEN
    3885                     IF ( ql(k,j,i) == 0.0_wp )  THEN
    3886                        k1 = 1.0_wp + 0.61_wp * q(k,j,i)
    3887                        k2 = 0.61_wp * pt(k,j,i)
    3888                     ELSE
    3889                        theta = pt(k,j,i) + d_exner(k) * lv_d_cp * ql(k,j,i)
    3890                        temp  = theta * exner(k)
    3891                        k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *                    &
    3892                                      ( q(k,j,i) - ql(k,j,i) ) *                &
    3893                             ( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) /          &
    3894                             ( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp *           &
    3895                             ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
    3896                        k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp )
    3897                    ENDIF
    3898                 ELSE IF ( cloud_droplets )  THEN
    3899                    k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
    3900                    k2 = 0.61_wp * pt(k,j,i)
    3901                 ENDIF
    3902 
    3903                 tend(k,j,i) = tend(k,j,i) + g /                                &
    3904                                 MERGE( vpt_reference, vpt(k,j,i),              &
    3905                                        use_single_reference_value ) *          &
    3906                                          ( k1 * surf_lsm_h%shf(m) +            &
    3907                                            k2 * surf_lsm_h%qsws(m)             &
    3908                                          ) * drho_air_zw(k-1)
    3909              ENDDO
    3910 !
    3911 !--          Treat horizontal urban surfaces
    3912              surf_s = surf_usm_h%start_index(j,i)
    3913              surf_e = surf_usm_h%end_index(j,i)
    3914              DO  m = surf_s, surf_e
    3915                 k = surf_usm_h%k(m)
    3916 
    3917                 IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets ) THEN
    3918                     k1 = 1.0_wp + 0.61_wp * q(k,j,i)
    3919                     k2 = 0.61_wp * pt(k,j,i)
    3920                 ELSE IF ( bulk_cloud_model )  THEN
    3921                     IF ( ql(k,j,i) == 0.0_wp )  THEN
    3922                        k1 = 1.0_wp + 0.61_wp * q(k,j,i)
    3923                        k2 = 0.61_wp * pt(k,j,i)
    3924                     ELSE
    3925                        theta = pt(k,j,i) + d_exner(k) * lv_d_cp * ql(k,j,i)
    3926                        temp  = theta * exner(k)
    3927                        k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *                    &
    3928                                      ( q(k,j,i) - ql(k,j,i) ) *                &
    3929                             ( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) /          &
    3930                             ( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp *           &
    3931                             ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
    3932                        k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp )
    3933                    ENDIF
    3934                 ELSE IF ( cloud_droplets )  THEN
    3935                    k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
    3936                    k2 = 0.61_wp * pt(k,j,i)
    3937                 ENDIF
    3938 
    3939                 tend(k,j,i) = tend(k,j,i) + g /                                &
    3940                                 MERGE( vpt_reference, vpt(k,j,i),              &
    3941                                        use_single_reference_value ) *          &
    3942                                          ( k1 * surf_usm_h%shf(m) +            &
    3943                                            k2 * surf_usm_h%qsws(m)             &
    3944                                          ) * drho_air_zw(k-1)
    3945              ENDDO
    3946 
    3947           ENDIF
    3948 
    3949           IF ( use_top_fluxes )  THEN
    3950              surf_s = surf_def_h(2)%start_index(j,i)
    3951              surf_e = surf_def_h(2)%end_index(j,i)
    3952              DO  m = surf_s, surf_e
    3953                 k = surf_def_h(2)%k(m)
    3954 
    3955 
    3956 
    3957                 IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets )  THEN
    39583728                   k1 = 1.0_wp + 0.61_wp * q(k,j,i)
    39593729                   k2 = 0.61_wp * pt(k,j,i)
     
    39653735                      theta = pt(k,j,i) + d_exner(k) * lv_d_cp * ql(k,j,i)
    39663736                      temp  = theta * exner(k)
    3967                       k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *                     &
    3968                                  ( q(k,j,i) - ql(k,j,i) ) *                    &
    3969                         ( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) /              &
    3970                         ( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp *               &
     3737                      k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *               &
     3738                                    ( q(k,j,i) - ql(k,j,i) ) *           &
     3739                           ( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) /     &
     3740                           ( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp *      &
     3741                           ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
     3742                      k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp )
     3743                   ENDIF
     3744                ELSE IF ( cloud_droplets )  THEN
     3745                   k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
     3746                   k2 = 0.61_wp * pt(k,j,i)
     3747                ENDIF
     3748
     3749                tmp_flux(k) = ( k1 * surf_lsm_h%shf(m) +                 &
     3750                                k2 * surf_lsm_h%qsws(m)                  &
     3751                              ) * drho_air_zw(k-1)
     3752             ENDDO
     3753!
     3754!--          Treat horizontal urban surfaces
     3755             surf_s = surf_usm_h%start_index(j,i)
     3756             surf_e = surf_usm_h%end_index(j,i)
     3757             DO  m = surf_s, surf_e
     3758                k = surf_usm_h%k(m)
     3759
     3760                IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets ) THEN
     3761                   k1 = 1.0_wp + 0.61_wp * q(k,j,i)
     3762                   k2 = 0.61_wp * pt(k,j,i)
     3763                ELSE IF ( bulk_cloud_model )  THEN
     3764                   IF ( ql(k,j,i) == 0.0_wp )  THEN
     3765                      k1 = 1.0_wp + 0.61_wp * q(k,j,i)
     3766                      k2 = 0.61_wp * pt(k,j,i)
     3767                   ELSE
     3768                      theta = pt(k,j,i) + d_exner(k) * lv_d_cp * ql(k,j,i)
     3769                      temp  = theta * exner(k)
     3770                      k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *               &
     3771                                    ( q(k,j,i) - ql(k,j,i) ) *           &
     3772                           ( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) /     &
     3773                           ( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp *      &
     3774                           ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
     3775                      k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp )
     3776                   ENDIF
     3777                ELSE IF ( cloud_droplets )  THEN
     3778                   k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
     3779                   k2 = 0.61_wp * pt(k,j,i)
     3780                ENDIF
     3781
     3782                tmp_flux(k) = ( k1 * surf_usm_h%shf(m) +                 &
     3783                                k2 * surf_usm_h%qsws(m)                  &
     3784                              ) * drho_air_zw(k-1)
     3785             ENDDO
     3786
     3787          ENDIF ! from IF ( use_surface_fluxes )  THEN
     3788
     3789          IF ( use_top_fluxes )  THEN
     3790
     3791             surf_s = surf_def_h(2)%start_index(j,i)
     3792             surf_e = surf_def_h(2)%end_index(j,i)
     3793             DO  m = surf_s, surf_e
     3794                k = surf_def_h(2)%k(m)
     3795
     3796                IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets ) THEN
     3797                   k1 = 1.0_wp + 0.61_wp * q(k,j,i)
     3798                   k2 = 0.61_wp * pt(k,j,i)
     3799                ELSE IF ( bulk_cloud_model )  THEN
     3800                   IF ( ql(k,j,i) == 0.0_wp )  THEN
     3801                      k1 = 1.0_wp + 0.61_wp * q(k,j,i)
     3802                      k2 = 0.61_wp * pt(k,j,i)
     3803                   ELSE
     3804                      theta = pt(k,j,i) + d_exner(k) * lv_d_cp * ql(k,j,i)
     3805                      temp  = theta * exner(k)
     3806                      k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *               &
     3807                                 ( q(k,j,i) - ql(k,j,i) ) *              &
     3808                        ( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) /        &
     3809                        ( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp *         &
    39713810                        ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
    39723811                      k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp )
     
    39773816                ENDIF
    39783817
    3979                 tend(k,j,i) = tend(k,j,i) + g /                                &
    3980                                 MERGE( vpt_reference, vpt(k,j,i),              &
    3981                                        use_single_reference_value ) *          &
    3982                             ( k1* surf_def_h(2)%shf(m) +                       &
    3983                               k2 * surf_def_h(2)%qsws(m)                       &
    3984                             ) * drho_air_zw(k)
     3818                tmp_flux(k) = ( k1 * surf_def_h(2)%shf(m) +              &
     3819                                k2 * surf_def_h(2)%qsws(m)               &
     3820                              ) * drho_air_zw(k)
     3821
     3822             ENDDO
     3823
     3824          ENDIF ! from IF ( use_top_fluxes )  THEN
     3825
     3826          IF ( .NOT. diss_production )  THEN
     3827
     3828!--          Compute temdency for TKE-production from shear
     3829             DO  k = nzb+1, nzt
     3830                flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_0(k,j,i),0) )
     3831                tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g / &
     3832                              MERGE( vpt_reference, vpt(k,j,i),          &
     3833                                     use_single_reference_value ) )
     3834             ENDDO
     3835
     3836          ELSE
     3837
     3838!--          RANS mode: Compute temdency for dissipation-rate-production from shear
     3839             DO  k = nzb+1, nzt
     3840                flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_0(k,j,i),0) )
     3841                tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g /   &
     3842                              MERGE( vpt_reference, vpt(k,j,i),          &
     3843                                     use_single_reference_value ) ) *    &
     3844                              diss(k,j,i)/( e(k,j,i) + 1.0E-20_wp ) *    &
     3845                              c_3
    39853846             ENDDO
    39863847
     
    39913852    ENDIF
    39923853
    3993   END SUBROUTINE production_e_ij
     3854 END SUBROUTINE production_e_ij
    39943855
    39953856
Note: See TracChangeset for help on using the changeset viewer.