- Timestamp:
- Oct 22, 2018 7:30:24 PM (6 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/turbulence_closure_mod.f90 (modified) (45 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/turbulence_closure_mod.f90
r3386 r3398 25 25 ! ----------------- 26 26 ! $Id$ 27 ! Refactored production_e and production_e_ij (removed code redundancy) 28 ! 29 ! 3386 2018-10-19 16:28:22Z gronemeier 27 30 ! Renamed tcm_prognostic to tcm_prognostic_equations 28 31 ! … … 227 230 REAL(wp) :: c_1 !< model constant for RANS mode 228 231 REAL(wp) :: c_2 !< model constant for RANS mode 229 !REAL(wp) :: c_3 !< model constant for RANS mode232 REAL(wp) :: c_3 !< model constant for RANS mode 230 233 REAL(wp) :: c_4 !< model constant for RANS mode 231 234 REAL(wp) :: l_max !< maximum length scale for Blackadar mixing length … … 234 237 235 238 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 closure239 (/ 0.55_wp, 1.44_wp, 1.92_wp, 1.44_wp, 0.0_wp /) !> default values fit for standard-tke-e closure 237 240 238 241 REAL(wp), DIMENSION(2) :: rans_const_sigma = & !< model constants for RANS mode, sigma values (sigma_e, sigma_diss) (namelist param) … … 423 426 c_1 = rans_const_c(1) 424 427 c_2 = rans_const_c(2) 425 !c_3 = rans_const_c(3) !> @todo clarify how to switch between different models428 c_3 = rans_const_c(3) !> @todo clarify how to switch between different models 426 429 c_4 = rans_const_c(4) 427 430 … … 2094 2097 IF ( rans_tke_e ) advec = tend 2095 2098 2096 CALL production_e 2099 CALL production_e( .FALSE. ) 2097 2100 2098 2101 ! … … 2547 2550 !> @todo Adjust production term in case of rans_tke_e simulation 2548 2551 !------------------------------------------------------------------------------! 2549 SUBROUTINE production_e 2552 SUBROUTINE production_e( diss_production ) 2550 2553 2551 2554 USE arrays_3d, & … … 2568 2571 2569 2572 IMPLICIT NONE 2573 2574 LOGICAL :: diss_production 2570 2575 2571 2576 INTEGER(iwp) :: i !< running index x-direction … … 2600 2605 REAL(wp), DIMENSION(nzb+1:nzt) :: dwdy !< Gradient of w-component in y-direction 2601 2606 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 2602 2608 2603 2609 … … 2824 2830 flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_0(k,j,i),flag_nr) ) 2825 2831 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 2827 2844 2828 2845 ENDDO … … 2844 2861 DO i = nxl, nxr 2845 2862 DO j = nys, nyn 2863 2846 2864 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) 2859 2866 ENDDO 2860 2867 ! … … 2867 2874 DO m = surf_s, surf_e 2868 2875 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) 2874 2877 ENDDO 2875 2878 ENDDO 2876 2877 2879 ENDIF 2878 2880 … … 2882 2884 DO m = surf_s, surf_e 2883 2885 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) 2889 2887 ENDDO 2890 2888 ENDIF 2891 2889 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 2892 2914 ENDDO 2893 2915 ENDDO … … 2899 2921 2900 2922 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) 2916 2924 ENDDO 2917 2925 … … 2924 2932 DO m = surf_s, surf_e 2925 2933 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) 2931 2935 ENDDO 2932 2936 ENDDO … … 2937 2941 DO m = surf_s, surf_e 2938 2942 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) 2944 2944 ENDDO 2945 2945 ! … … 2949 2949 DO m = surf_s, surf_e 2950 2950 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) 2956 2952 ENDDO 2957 2953 ENDIF … … 2962 2958 DO m = surf_s, surf_e 2963 2959 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) 2969 2961 ENDDO 2970 2962 ENDIF 2971 2963 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 2972 2988 ENDDO 2973 2989 ENDDO … … 2981 2997 2982 2998 DO k = nzb+1, nzt 2983 ! 2984 !-- Flag 9 is used to mask top fluxes, flag 30 to mask 2985 !-- surface fluxes 2999 2986 3000 IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets ) THEN 2987 3001 k1 = 1.0_wp + 0.61_wp * q(k,j,i) 2988 3002 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) 3002 3007 ELSE IF ( bulk_cloud_model ) THEN 3003 3008 IF ( ql(k,j,i) == 0.0_wp ) THEN … … 3014 3019 k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp ) 3015 3020 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) 3029 3025 ELSE IF ( cloud_droplets ) THEN 3030 3026 k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i) 3031 3027 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) 3046 3033 ENDIF 3047 3034 … … 3079 3066 ENDIF 3080 3067 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) 3087 3071 ENDDO 3088 3072 ENDDO … … 3116 3100 ENDIF 3117 3101 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) 3124 3105 ENDDO 3125 3106 ! … … 3152 3133 ENDIF 3153 3134 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) 3160 3138 ENDDO 3161 3139 … … 3191 3169 ENDIF 3192 3170 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) 3199 3174 3200 3175 ENDDO 3201 3176 3202 3177 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 3203 3202 3204 3203 ENDDO … … 3227 3226 3228 3227 USE control_parameters, & 3229 ONLY: cloud_droplets, constant_flux_layer, neutral, &3228 ONLY: cloud_droplets, constant_flux_layer, neutral, & 3230 3229 rho_reference, use_single_reference_value, use_surface_fluxes, & 3231 3230 use_top_fluxes … … 3239 3238 USE surface_mod, & 3240 3239 ONLY : surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, & 3241 surf_usm_v 3240 surf_usm_v 3242 3241 3243 3242 IMPLICIT NONE … … 3252 3251 INTEGER(iwp) :: surf_e !< end index of surface elements at given i-j position 3253 3252 INTEGER(iwp) :: surf_s !< start index of surface elements at given i-j position 3253 INTEGER(iwp) :: flag_nr !< number of masking flag 3254 3254 3255 3255 REAL(wp) :: def !< … … 3260 3260 REAL(wp) :: theta !< 3261 3261 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 3263 3263 REAL(wp) :: usvs !< momentum flux u"v" 3264 3264 REAL(wp) :: vsus !< momentum flux v"u" … … 3266 3266 REAL(wp) :: wsvs !< momentum flux w"v" 3267 3267 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 3279 3308 3280 3309 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' 3306 3316 ! 3307 3317 !-- 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 3309 3320 DO l = 0, 1 3310 !3311 !-- Default surfaces3312 3321 surf_s = surf_def_v(l)%start_index(j,i) 3313 3322 surf_e = surf_def_v(l)%end_index(j,i) … … 3317 3326 wsvs = surf_def_v(l)%mom_flux_tke(1,m) 3318 3327 3319 km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp 3328 km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp & 3320 3329 * 0.5_wp * dy 3321 3330 ! 3322 3331 !-- -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 ) ) 3325 3334 dudy(k) = sign_dir * usvs / ( km_neutral + 1E-10_wp ) 3326 3335 dwdy(k) = sign_dir * wsvs / ( km_neutral + 1E-10_wp ) 3327 ENDDO 3336 ENDDO 3328 3337 ! 3329 3338 !-- Natural surfaces … … 3335 3344 wsvs = surf_lsm_v(l)%mom_flux_tke(1,m) 3336 3345 3337 km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp 3346 km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp & 3338 3347 * 0.5_wp * dy 3339 3348 ! 3340 3349 !-- -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 ) ) 3343 3352 dudy(k) = sign_dir * usvs / ( km_neutral + 1E-10_wp ) 3344 3353 dwdy(k) = sign_dir * wsvs / ( km_neutral + 1E-10_wp ) 3345 ENDDO 3354 ENDDO 3346 3355 ! 3347 3356 !-- Urban surfaces … … 3353 3362 wsvs = surf_usm_v(l)%mom_flux_tke(1,m) 3354 3363 3355 km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp 3364 km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp & 3356 3365 * 0.5_wp * dy 3357 3366 ! 3358 3367 !-- -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 ) ) 3361 3370 dudy(k) = sign_dir * usvs / ( km_neutral + 1E-10_wp ) 3362 3371 dwdy(k) = sign_dir * wsvs / ( km_neutral + 1E-10_wp ) 3363 ENDDO 3372 ENDDO 3364 3373 ENDDO 3365 3374 ! 3366 !-- Compute gradients at east- and west-facing walls3375 !-- Compute gradients at east- and west-facing walls 3367 3376 DO l = 2, 3 3368 !3369 !-- Default surfaces3370 3377 surf_s = surf_def_v(l)%start_index(j,i) 3371 3378 surf_e = surf_def_v(l)%end_index(j,i) 3372 3379 DO m = surf_s, surf_e 3373 k 3374 vsus 3375 wsus 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 & 3378 3385 * 0.5_wp * dx 3379 3386 ! 3380 3387 !-- -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 ) ) 3383 3390 dvdx(k) = sign_dir * vsus / ( km_neutral + 1E-10_wp ) 3384 3391 dwdx(k) = sign_dir * wsus / ( km_neutral + 1E-10_wp ) 3385 ENDDO 3386 ! 3387 !-- Natural surfaces3392 ENDDO 3393 ! 3394 !-- Natural surfaces 3388 3395 surf_s = surf_lsm_v(l)%start_index(j,i) 3389 3396 surf_e = surf_lsm_v(l)%end_index(j,i) 3390 3397 DO m = surf_s, surf_e 3391 k 3392 vsus 3393 wsus 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 & 3396 3403 * 0.5_wp * dx 3397 3404 ! 3398 3405 !-- -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 ) ) 3401 3408 dvdx(k) = sign_dir * vsus / ( km_neutral + 1E-10_wp ) 3402 3409 dwdx(k) = sign_dir * wsus / ( km_neutral + 1E-10_wp ) 3403 ENDDO 3404 ! 3405 !-- Urban surfaces3410 ENDDO 3411 ! 3412 !-- Urban surfaces 3406 3413 surf_s = surf_usm_v(l)%start_index(j,i) 3407 3414 surf_e = surf_usm_v(l)%end_index(j,i) 3408 3415 DO m = surf_s, surf_e 3409 k 3410 vsus 3411 wsus 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 & 3414 3421 * 0.5_wp * dx 3415 3422 ! 3416 3423 !-- -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 ) ) 3419 3426 dvdx(k) = sign_dir * vsus / ( km_neutral + 1E-10_wp ) 3420 3427 dwdx(k) = sign_dir * wsus / ( km_neutral + 1E-10_wp ) 3421 ENDDO 3428 ENDDO 3422 3429 ENDDO 3423 3430 ! 3424 !-- Compute gradients at upward-facing walls, first for 3425 !-- non-natural default surfaces 3431 !-- Compute gradients at upward-facing surfaces 3426 3432 surf_s = surf_def_h(0)%start_index(j,i) 3427 3433 surf_e = surf_def_h(0)%end_index(j,i) … … 3429 3435 k = surf_def_h(0)%k(m) 3430 3436 ! 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 3433 3439 !-- 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 3435 3441 !-- 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) 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) 3438 3444 3439 3445 ENDDO … … 3445 3451 k = surf_lsm_h%k(m) 3446 3452 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 3449 3456 ENDDO 3450 3457 ! … … 3455 3462 k = surf_usm_h%k(m) 3456 3463 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 3459 3467 ENDDO 3460 3468 ! 3461 !-- Compute gradients at downward-facing walls, only for 3469 !-- Compute gradients at downward-facing walls, only for 3462 3470 !-- non-natural default surfaces 3463 3471 surf_s = surf_def_h(1)%start_index(j,i) … … 3466 3474 k = surf_def_h(1)%k(m) 3467 3475 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) 3470 3478 3471 3479 ENDDO 3472 3480 3473 ! IF ( .NOT. rans_tke_e ) THEN3474 3475 DO k = nzb+1, nzt3476 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_wp3485 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 * flag3489 3490 ENDDO3491 3492 ! ELSE3493 !3494 ! DO k = nzb+1, nzt3495 ! !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_wp3510 !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 * flag3514 !3515 ! ENDDO3516 !3517 ! ENDIF3518 3519 ELSE ! not constant_flux_layer3520 3521 ! IF ( .NOT. rans_tke_e ) THEN3522 !3523 !-- Calculate TKE production by shear. Here, no additional3524 !-- wall-bounded code is considered.3525 !-- Why?3526 DO k = nzb+1, nzt3527 3528 dudx(k) = ( u(k,j,i+1) - u(k,j,i) ) * ddx3529 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) ) * ddy3531 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) ) * ddx3536 dvdy(k) = ( v(k,j+1,i) - v(k,j,i) ) * ddy3537 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) ) * ddx3542 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) ) * ddy3544 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_wp3554 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 * flag3558 3559 ENDDO3560 3561 ! ELSE3562 !3563 ! DO k = nzb+1, nzt3564 !3565 ! dudx(k) = ( u(k,j,i+1) - u(k,j,i) ) * ddx3566 ! 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) ) * ddy3568 ! 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) ) * ddx3573 ! dvdy(k) = ( v(k,j+1,i) - v(k,j,i) ) * ddy3574 ! 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) ) * ddx3579 ! 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) ) * ddy3581 ! 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_wp3597 !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 * flag3601 !3602 ! ENDDO3603 !3604 ! ENDIF3605 3606 3481 ENDIF 3607 3482 3608 IF ( .NOT. diss_production ) THEN3609 ! 3610 !-- Production term in case of TKE production 3611 DO k = nzb+1, nzt3612 tend(k,j,i) = tend(k,j,i) + tend_temp(k)3613 ENDDO3614 ELSE3615 ! 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 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 3635 3510 3636 3511 ! … … 3642 3517 IF ( ocean_mode ) THEN 3643 3518 ! 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 3646 3522 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) 3660 3577 ENDDO 3661 3578 … … 3668 3585 DO m = surf_s, surf_e 3669 3586 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) 3675 3588 ENDDO 3676 3589 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 3678 3606 ENDIF 3679 3607 … … 3683 3611 DO m = surf_s, surf_e 3684 3612 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) 3690 3614 ENDDO 3691 3615 ENDIF 3692 3616 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 ) ) 3727 3625 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_e3733 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 3739 3637 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 3752 3639 ENDIF 3753 3640 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 3770 3644 3771 3645 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 3775 3648 k1 = 1.0_wp + 0.61_wp * q(k,j,i) 3776 3649 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) 3790 3654 ELSE IF ( bulk_cloud_model ) THEN 3791 3655 IF ( ql(k,j,i) == 0.0_wp ) THEN … … 3795 3659 theta = pt(k,j,i) + d_exner(k) * lv_d_cp * ql(k,j,i) 3796 3660 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 * & 3801 3665 ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) ) 3802 3666 k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp ) 3803 3667 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) 3816 3672 ELSE IF ( cloud_droplets ) THEN 3817 3673 k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i) 3818 3674 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) 3832 3680 ENDIF 3681 3833 3682 ENDDO 3834 3683 3835 3684 IF ( use_surface_fluxes ) THEN 3836 3685 ! 3837 !-- Treat horizontal default surfaces , up- and downward-facing3686 !-- Treat horizontal default surfaces 3838 3687 DO l = 0, 1 3839 3688 surf_s = surf_def_h(l)%start_index(j,i) … … 3850 3699 k2 = 0.61_wp * pt(k,j,i) 3851 3700 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 ) 3860 3709 ENDIF 3861 3710 ELSE IF ( cloud_droplets ) THEN … … 3864 3713 ENDIF 3865 3714 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) 3872 3718 ENDDO 3873 3719 ENDDO … … 3880 3726 3881 3727 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 ) THEN3885 IF ( ql(k,j,i) == 0.0_wp ) THEN3886 k1 = 1.0_wp + 0.61_wp * q(k,j,i)3887 k2 = 0.61_wp * pt(k,j,i)3888 ELSE3889 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 ENDIF3898 ELSE IF ( cloud_droplets ) THEN3899 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 ENDIF3902 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 ENDDO3910 !3911 !-- Treat horizontal urban surfaces3912 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_e3915 k = surf_usm_h%k(m)3916 3917 IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets ) THEN3918 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 ) THEN3921 IF ( ql(k,j,i) == 0.0_wp ) THEN3922 k1 = 1.0_wp + 0.61_wp * q(k,j,i)3923 k2 = 0.61_wp * pt(k,j,i)3924 ELSE3925 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 ENDIF3934 ELSE IF ( cloud_droplets ) THEN3935 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 ENDIF3938 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 ENDDO3946 3947 ENDIF3948 3949 IF ( use_top_fluxes ) THEN3950 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_e3953 k = surf_def_h(2)%k(m)3954 3955 3956 3957 IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets ) THEN3958 3728 k1 = 1.0_wp + 0.61_wp * q(k,j,i) 3959 3729 k2 = 0.61_wp * pt(k,j,i) … … 3965 3735 theta = pt(k,j,i) + d_exner(k) * lv_d_cp * ql(k,j,i) 3966 3736 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 * & 3971 3810 ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) ) 3972 3811 k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp ) … … 3977 3816 ENDIF 3978 3817 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 3985 3846 ENDDO 3986 3847 … … 3991 3852 ENDIF 3992 3853 3993 3854 END SUBROUTINE production_e_ij 3994 3855 3995 3856
Note: See TracChangeset
for help on using the changeset viewer.
