!> @file production_e.f90 !------------------------------------------------------------------------------! ! This file is part of PALM. ! ! PALM is free software: you can redistribute it and/or modify it under the ! terms of the GNU General Public License as published by the Free Software ! Foundation, either version 3 of the License, or (at your option) any later ! version. ! ! PALM is distributed in the hope that it will be useful, but WITHOUT ANY ! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR ! A PARTICULAR PURPOSE. See the GNU General Public License for more details. ! ! You should have received a copy of the GNU General Public License along with ! PALM. If not, see . ! ! Copyright 1997-2017 Leibniz Universitaet Hannover !------------------------------------------------------------------------------! ! ! Current revisions: ! ----------------- ! ! ! Former revisions: ! ----------------- ! $Id: production_e.f90 2233 2017-05-30 18:08:54Z scharf $ ! ! 2232 2017-05-30 17:47:52Z suehring ! Adjustments to new surface concept ! ! 2126 2017-01-20 15:54:21Z raasch ! density in ocean case replaced by potential density ! ! 2118 2017-01-17 16:38:49Z raasch ! OpenACC version of subroutine removed ! ! 2031 2016-10-21 15:11:58Z knoop ! renamed variable rho to rho_ocean ! ! 2000 2016-08-20 18:09:15Z knoop ! Forced header and separation lines into 80 columns ! ! 1873 2016-04-18 14:50:06Z maronga ! Module renamed (removed _mod) ! ! ! 1850 2016-04-08 13:29:27Z maronga ! Module renamed ! ! ! 1691 2015-10-26 16:17:44Z maronga ! Renamed prandtl_layer to constant_flux_layer. ! ! 1682 2015-10-07 23:56:08Z knoop ! Code annotations made doxygen readable ! ! 1374 2014-04-25 12:55:07Z raasch ! nzb_s_outer removed from acc-present-list ! ! 1353 2014-04-08 15:21:23Z heinze ! REAL constants provided with KIND-attribute ! ! 1342 2014-03-26 17:04:47Z kanani ! REAL constants defined as wp-kind ! ! 1320 2014-03-20 08:40:49Z raasch ! ONLY-attribute added to USE-statements, ! kind-parameters added to all INTEGER and REAL declaration statements, ! kinds are defined in new module kinds, ! old module precision_kind is removed, ! revision history before 2012 removed, ! comment fields (!:) to be used for variable explanations added to ! all variable declaration statements ! ! 1257 2013-11-08 15:18:40Z raasch ! openacc loop and loop vector clauses removed, declare create moved after ! the FORTRAN declaration statement ! ! 1179 2013-06-14 05:57:58Z raasch ! use_reference renamed use_single_reference_value ! ! 1128 2013-04-12 06:19:32Z raasch ! loop index bounds in accelerator version replaced by i_left, i_right, j_south, ! j_north ! ! 1036 2012-10-22 13:43:42Z raasch ! code put under GPL (PALM 3.9) ! ! 1015 2012-09-27 09:23:24Z raasch ! accelerator version (*_acc) added ! ! 1007 2012-09-19 14:30:36Z franke ! Bugfix: calculation of buoyancy production has to consider the liquid water ! mixing ratio in case of cloud droplets ! ! 940 2012-07-09 14:31:00Z raasch ! TKE production by buoyancy can be switched off in case of runs with pure ! neutral stratification ! ! Revision 1.1 1997/09/19 07:45:35 raasch ! Initial revision ! ! ! Description: ! ------------ !> Production terms (shear + buoyancy) of the TKE. !> @warning The case with constant_flux_layer = F and use_surface_fluxes = T is !> not considered well! !------------------------------------------------------------------------------! MODULE production_e_mod USE kinds PRIVATE PUBLIC production_e, production_e_init INTERFACE production_e MODULE PROCEDURE production_e MODULE PROCEDURE production_e_ij END INTERFACE production_e INTERFACE production_e_init MODULE PROCEDURE production_e_init END INTERFACE production_e_init CONTAINS !------------------------------------------------------------------------------! ! Description: ! ------------ !> Call for all grid points !------------------------------------------------------------------------------! SUBROUTINE production_e USE arrays_3d, & ONLY: ddzw, dd2zu, kh, km, prho, pt, q, ql, tend, u, v, vpt, w USE cloud_parameters, & ONLY: l_d_cp, l_d_r, pt_d_t, t_d_pt USE control_parameters, & ONLY: cloud_droplets, cloud_physics, constant_flux_layer, g, & humidity, kappa, neutral, ocean, pt_reference, & rho_reference, use_single_reference_value, & use_surface_fluxes, use_top_fluxes USE grid_variables, & ONLY: ddx, dx, ddy, dy USE indices, & ONLY: nxl, nxr, nys, nyn, nzb, nzt, wall_flags_0 USE surface_mod, & ONLY : surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, & surf_usm_v IMPLICIT NONE INTEGER(iwp) :: i !< running index x-direction INTEGER(iwp) :: j !< running index y-direction INTEGER(iwp) :: k !< running index z-direction INTEGER(iwp) :: l !< running index for different surface type orientation INTEGER(iwp) :: m !< running index surface elements INTEGER(iwp) :: surf_e !< end index of surface elements at given i-j position INTEGER(iwp) :: surf_s !< start index of surface elements at given i-j position REAL(wp) :: def !< REAL(wp) :: flag !< flag to mask topography REAL(wp) :: k1 !< REAL(wp) :: k2 !< REAL(wp) :: km_neutral !< diffusion coefficient assuming neutral conditions - used to compute shear production at surfaces REAL(wp) :: theta !< REAL(wp) :: temp !< REAL(wp) :: sign_dir !< sign of wall-tke flux, depending on wall orientation REAL(wp) :: usvs !< momentum flux u"v" REAL(wp) :: vsus !< momentum flux v"u" REAL(wp) :: wsus !< momentum flux w"u" REAL(wp) :: wsvs !< momentum flux w"v" REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: dudx !< Gradient of u-component in x-direction REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: dudy !< Gradient of u-component in y-direction REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: dudz !< Gradient of u-component in z-direction REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: dvdx !< Gradient of v-component in x-direction REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: dvdy !< Gradient of v-component in y-direction REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: dvdz !< Gradient of v-component in z-direction REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: dwdx !< Gradient of w-component in x-direction REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: dwdy !< Gradient of w-component in y-direction REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: dwdz !< Gradient of w-component in z-direction DO i = nxl, nxr IF ( constant_flux_layer ) THEN ! !-- Calculate TKE production by shear. Calculate gradients at all grid !-- points first, gradients at surface-bounded grid points will be !-- overwritten further below. DO j = nys, nyn DO k = nzb+1, nzt dudx(k,j) = ( u(k,j,i+1) - u(k,j,i) ) * ddx dudy(k,j) = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - & u(k,j-1,i) - u(k,j-1,i+1) ) * ddy dudz(k,j) = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - & u(k-1,j,i) - u(k-1,j,i+1) ) * & dd2zu(k) dvdx(k,j) = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - & v(k,j,i-1) - v(k,j+1,i-1) ) * ddx dvdy(k,j) = ( v(k,j+1,i) - v(k,j,i) ) * ddy dvdz(k,j) = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - & v(k-1,j,i) - v(k-1,j+1,i) ) * & dd2zu(k) dwdx(k,j) = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - & w(k,j,i-1) - w(k-1,j,i-1) ) * ddx dwdy(k,j) = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - & w(k,j-1,i) - w(k-1,j-1,i) ) * ddy dwdz(k,j) = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) ENDDO ENDDO ! !-- Position beneath wall !-- (2) - Will allways be executed. !-- 'bottom and wall: use u_0,v_0 and wall functions' DO j = nys, nyn ! !-- Compute gradients at north- and south-facing surfaces. !-- First, for default surfaces, then for urban surfaces. !-- Note, so far no natural vertical surfaces implemented DO l = 0, 1 surf_s = surf_def_v(l)%start_index(j,i) surf_e = surf_def_v(l)%end_index(j,i) DO m = surf_s, surf_e k = surf_def_v(l)%k(m) usvs = surf_def_v(l)%mom_flux_tke(0,m) wsvs = surf_def_v(l)%mom_flux_tke(1,m) km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp & * 0.5_wp * dy ! !-- -1.0 for right-facing wall, 1.0 for left-facing wall sign_dir = MERGE( 1.0_wp, -1.0_wp, & BTEST( wall_flags_0(k,j-1,i), 0 ) ) dudy(k,j) = sign_dir * usvs / ( km_neutral + 1E-10_wp ) dwdy(k,j) = sign_dir * wsvs / ( km_neutral + 1E-10_wp ) ENDDO ! !-- Natural surfaces surf_s = surf_lsm_v(l)%start_index(j,i) surf_e = surf_lsm_v(l)%end_index(j,i) DO m = surf_s, surf_e k = surf_lsm_v(l)%k(m) usvs = surf_lsm_v(l)%mom_flux_tke(0,m) wsvs = surf_lsm_v(l)%mom_flux_tke(1,m) km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp & * 0.5_wp * dy ! !-- -1.0 for right-facing wall, 1.0 for left-facing wall sign_dir = MERGE( 1.0_wp, -1.0_wp, & BTEST( wall_flags_0(k,j-1,i), 0 ) ) dudy(k,j) = sign_dir * usvs / ( km_neutral + 1E-10_wp ) dwdy(k,j) = sign_dir * wsvs / ( km_neutral + 1E-10_wp ) ENDDO ! !-- Urban surfaces surf_s = surf_usm_v(l)%start_index(j,i) surf_e = surf_usm_v(l)%end_index(j,i) DO m = surf_s, surf_e k = surf_usm_v(l)%k(m) usvs = surf_usm_v(l)%mom_flux_tke(0,m) wsvs = surf_usm_v(l)%mom_flux_tke(1,m) km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp & * 0.5_wp * dy ! !-- -1.0 for right-facing wall, 1.0 for left-facing wall sign_dir = MERGE( 1.0_wp, -1.0_wp, & BTEST( wall_flags_0(k,j-1,i), 0 ) ) dudy(k,j) = sign_dir * usvs / ( km_neutral + 1E-10_wp ) dwdy(k,j) = sign_dir * wsvs / ( km_neutral + 1E-10_wp ) ENDDO ENDDO ! !-- Compute gradients at east- and west-facing walls DO l = 2, 3 surf_s = surf_def_v(l)%start_index(j,i) surf_e = surf_def_v(l)%end_index(j,i) DO m = surf_s, surf_e k = surf_def_v(l)%k(m) vsus = surf_def_v(l)%mom_flux_tke(0,m) wsus = surf_def_v(l)%mom_flux_tke(1,m) km_neutral = kappa * ( vsus**2 + wsus**2 )**0.25_wp & * 0.5_wp * dx ! !-- -1.0 for right-facing wall, 1.0 for left-facing wall sign_dir = MERGE( 1.0_wp, -1.0_wp, & BTEST( wall_flags_0(k,j,i-1), 0 ) ) dvdx(k,j) = sign_dir * vsus / ( km_neutral + 1E-10_wp ) dwdx(k,j) = sign_dir * wsus / ( km_neutral + 1E-10_wp ) ENDDO ! !-- Natural surfaces surf_s = surf_lsm_v(l)%start_index(j,i) surf_e = surf_lsm_v(l)%end_index(j,i) DO m = surf_s, surf_e k = surf_lsm_v(l)%k(m) vsus = surf_lsm_v(l)%mom_flux_tke(0,m) wsus = surf_lsm_v(l)%mom_flux_tke(1,m) km_neutral = kappa * ( vsus**2 + wsus**2 )**0.25_wp & * 0.5_wp * dx ! !-- -1.0 for right-facing wall, 1.0 for left-facing wall sign_dir = MERGE( 1.0_wp, -1.0_wp, & BTEST( wall_flags_0(k,j,i-1), 0 ) ) dvdx(k,j) = sign_dir * vsus / ( km_neutral + 1E-10_wp ) dwdx(k,j) = sign_dir * wsus / ( km_neutral + 1E-10_wp ) ENDDO ! !-- Urban surfaces surf_s = surf_usm_v(l)%start_index(j,i) surf_e = surf_usm_v(l)%end_index(j,i) DO m = surf_s, surf_e k = surf_usm_v(l)%k(m) vsus = surf_usm_v(l)%mom_flux_tke(0,m) wsus = surf_usm_v(l)%mom_flux_tke(1,m) km_neutral = kappa * ( vsus**2 + wsus**2 )**0.25_wp & * 0.5_wp * dx ! !-- -1.0 for right-facing wall, 1.0 for left-facing wall sign_dir = MERGE( 1.0_wp, -1.0_wp, & BTEST( wall_flags_0(k,j,i-1), 0 ) ) dvdx(k,j) = sign_dir * vsus / ( km_neutral + 1E-10_wp ) dwdx(k,j) = sign_dir * wsus / ( km_neutral + 1E-10_wp ) ENDDO ENDDO ! !-- Compute gradients at upward-facing surfaces surf_s = surf_def_h(0)%start_index(j,i) surf_e = surf_def_h(0)%end_index(j,i) DO m = surf_s, surf_e k = surf_def_h(0)%k(m) ! !-- Please note, actually, an interpolation of u_0 and v_0 !-- onto the grid center would be required. However, this !-- would require several data transfers between 2D-grid and !-- wall type. The effect of this missing interpolation is !-- negligible. (See also production_e_init). dudz(k,j) = ( u(k+1,j,i) - surf_def_h(0)%u_0(m) ) * dd2zu(k) dvdz(k,j) = ( v(k+1,j,i) - surf_def_h(0)%v_0(m) ) * dd2zu(k) ENDDO ! !-- Natural surfaces surf_s = surf_lsm_h%start_index(j,i) surf_e = surf_lsm_h%end_index(j,i) DO m = surf_s, surf_e k = surf_lsm_h%k(m) ! !-- Please note, actually, an interpolation of u_0 and v_0 !-- onto the grid center would be required. However, this !-- would require several data transfers between 2D-grid and !-- wall type. The effect of this missing interpolation is !-- negligible. (See also production_e_init). dudz(k,j) = ( u(k+1,j,i) - surf_lsm_h%u_0(m) ) * dd2zu(k) dvdz(k,j) = ( v(k+1,j,i) - surf_lsm_h%v_0(m) ) * dd2zu(k) ENDDO ! !-- Urban surfaces surf_s = surf_usm_h%start_index(j,i) surf_e = surf_usm_h%end_index(j,i) DO m = surf_s, surf_e k = surf_usm_h%k(m) ! !-- Please note, actually, an interpolation of u_0 and v_0 !-- onto the grid center would be required. However, this !-- would require several data transfers between 2D-grid and !-- wall type. The effect of this missing interpolation is !-- negligible. (See also production_e_init). dudz(k,j) = ( u(k+1,j,i) - surf_usm_h%u_0(m) ) * dd2zu(k) dvdz(k,j) = ( v(k+1,j,i) - surf_usm_h%v_0(m) ) * dd2zu(k) ENDDO ! !-- Compute gradients at downward-facing walls, only for !-- non-natural default surfaces surf_s = surf_def_h(1)%start_index(j,i) surf_e = surf_def_h(1)%end_index(j,i) DO m = surf_s, surf_e k = surf_def_h(1)%k(m) ! !-- Please note, actually, an interpolation of u_0 and v_0 !-- onto the grid center would be required. However, this !-- would require several data transfers between 2D-grid and !-- wall type. The effect of this missing interpolation is !-- negligible. (See also production_e_init). dudz(k,j) = ( surf_def_h(1)%u_0(m) - u(k-1,j,i) ) * dd2zu(k) dvdz(k,j) = ( surf_def_h(1)%v_0(m) - v(k-1,j,i) ) * dd2zu(k) ENDDO ENDDO DO j = nys, nyn DO k = nzb+1, nzt def = 2.0_wp * ( dudx(k,j)**2 + dvdy(k,j)**2 + dwdz(k,j)**2 ) + & dudy(k,j)**2 + dvdx(k,j)**2 + dwdx(k,j)**2 + & dwdy(k,j)**2 + dudz(k,j)**2 + dvdz(k,j)**2 + & 2.0_wp * ( dvdx(k,j)*dudy(k,j) + dwdx(k,j)*dudz(k,j) + & dwdy(k,j)*dvdz(k,j) ) IF ( def < 0.0_wp ) def = 0.0_wp flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def * flag ENDDO ENDDO ELSEIF ( use_surface_fluxes ) THEN DO j = nys, nyn ! !-- Calculate TKE production by shear. Here, no additional !-- wall-bounded code is considered. !-- Why? DO k = nzb+1, nzt dudx(k,j) = ( u(k,j,i+1) - u(k,j,i) ) * ddx dudy(k,j) = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - & u(k,j-1,i) - u(k,j-1,i+1) ) * ddy dudz(k,j) = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - & u(k-1,j,i) - u(k-1,j,i+1) ) * & dd2zu(k) dvdx(k,j) = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - & v(k,j,i-1) - v(k,j+1,i-1) ) * ddx dvdy(k,j) = ( v(k,j+1,i) - v(k,j,i) ) * ddy dvdz(k,j) = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - & v(k-1,j,i) - v(k-1,j+1,i) ) * & dd2zu(k) dwdx(k,j) = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - & w(k,j,i-1) - w(k-1,j,i-1) ) * ddx dwdy(k,j) = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - & w(k,j-1,i) - w(k-1,j-1,i) ) * ddy dwdz(k,j) = ( w(k,j,i) - w(k-1,j,i) ) * & ddzw(k) def = 2.0_wp * ( & dudx(k,j)**2 + dvdy(k,j)**2 + dwdz(k,j)**2 & ) + & dudy(k,j)**2 + dvdx(k,j)**2 + dwdx(k,j)**2 + & dwdy(k,j)**2 + dudz(k,j)**2 + dvdz(k,j)**2 + & 2.0_wp * ( & dvdx(k,j)*dudy(k,j) + dwdx(k,j)*dudz(k,j) + & dwdy(k,j)*dvdz(k,j) & ) IF ( def < 0.0_wp ) def = 0.0_wp flag = MERGE( 1.0_wp, 0.0_wp, & BTEST( wall_flags_0(k,j,i), 29 ) ) tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def * flag ENDDO ENDDO ENDIF ! !-- If required, calculate TKE production by buoyancy IF ( .NOT. neutral ) THEN IF ( .NOT. humidity ) THEN IF ( use_single_reference_value ) THEN IF ( ocean ) THEN ! !-- So far in the ocean no special treatment of density flux !-- in the bottom and top surface layer DO j = nys, nyn DO k = nzb+1, nzt tend(k,j,i) = tend(k,j,i) + & kh(k,j,i) * g / rho_reference * & ( prho(k+1,j,i) - prho(k-1,j,i) ) * & dd2zu(k) * & MERGE( 1.0_wp, 0.0_wp, & BTEST( wall_flags_0(k,j,i), 0 ) & ) ENDDO ENDDO ELSE DO j = nys, nyn DO k = nzb+1, nzt ! !-- Flag 9 is used to mask top fluxes, flag 30 to mask !-- surface fluxes tend(k,j,i) = tend(k,j,i) - & kh(k,j,i) * g / pt_reference * & ( pt(k+1,j,i) - pt(k-1,j,i) ) * & dd2zu(k) * & MERGE( 1.0_wp, 0.0_wp, & BTEST( wall_flags_0(k,j,i), 30 ) & ) * & MERGE( 1.0_wp, 0.0_wp, & BTEST( wall_flags_0(k,j,i), 9 ) & ) ENDDO IF ( use_surface_fluxes ) THEN ! !-- Default surfaces, up- and downward-facing DO l = 0, 1 surf_s = surf_def_h(l)%start_index(j,i) surf_e = surf_def_h(l)%end_index(j,i) DO m = surf_s, surf_e k = surf_def_h(l)%k(m) tend(k,j,i) = tend(k,j,i) + g / pt_reference & * surf_def_h(l)%shf(m) ENDDO ENDDO ! !-- Natural surfaces surf_s = surf_lsm_h%start_index(j,i) surf_e = surf_lsm_h%end_index(j,i) DO m = surf_s, surf_e k = surf_lsm_h%k(m) tend(k,j,i) = tend(k,j,i) + g / pt_reference & * surf_lsm_h%shf(m) ENDDO ! !-- Urban surfaces surf_s = surf_usm_h%start_index(j,i) surf_e = surf_usm_h%end_index(j,i) DO m = surf_s, surf_e k = surf_usm_h%k(m) tend(k,j,i) = tend(k,j,i) + g / pt_reference & * surf_usm_h%shf(m) ENDDO ENDIF IF ( use_top_fluxes ) THEN surf_s = surf_def_h(2)%start_index(j,i) surf_e = surf_def_h(2)%end_index(j,i) DO m = surf_s, surf_e k = surf_def_h(2)%k(m) tend(k,j,i) = tend(k,j,i) + g / pt_reference * & surf_def_h(2)%shf(m) ENDDO ENDIF ENDDO ENDIF ELSE IF ( ocean ) THEN ! !-- So far in the ocean no special treatment of density flux !-- in the bottom and top surface layer DO j = nys, nyn DO k = nzb+1, nzt tend(k,j,i) = tend(k,j,i) + & kh(k,j,i) * g / prho(k,j,i) * & ( prho(k+1,j,i) - prho(k-1,j,i) ) * & dd2zu(k) * & MERGE( 1.0_wp, 0.0_wp, & BTEST( wall_flags_0(k,j,i), 0 ) & ) ENDDO ENDDO ELSE DO j = nys, nyn DO k = nzb+1, nzt ! !-- Flag 9 is used to mask top fluxes, flag 30 to mask !-- surface fluxes tend(k,j,i) = tend(k,j,i) - & kh(k,j,i) * g / pt(k,j,i) * & ( pt(k+1,j,i) - pt(k-1,j,i) ) * & dd2zu(k) * & MERGE( 1.0_wp, 0.0_wp, & BTEST( wall_flags_0(k,j,i), 30 ) & ) * & MERGE( 1.0_wp, 0.0_wp, & BTEST( wall_flags_0(k,j,i), 9 ) & ) ENDDO IF ( use_surface_fluxes ) THEN ! !-- Default surfaces, up- and downwrd-facing DO l = 0, 1 surf_s = surf_def_h(l)%start_index(j,i) surf_e = surf_def_h(l)%end_index(j,i) DO m = surf_s, surf_e k = surf_def_h(l)%k(m) tend(k,j,i) = tend(k,j,i) + g / pt_reference & * surf_def_h(l)%shf(m) ENDDO ENDDO ! !-- Natural surfaces surf_s = surf_lsm_h%start_index(j,i) surf_e = surf_lsm_h%end_index(j,i) DO m = surf_s, surf_e k = surf_lsm_h%k(m) tend(k,j,i) = tend(k,j,i) + g / pt_reference & * surf_lsm_h%shf(m) ENDDO ! !-- Urban surfaces surf_s = surf_usm_h%start_index(j,i) surf_e = surf_usm_h%end_index(j,i) DO m = surf_s, surf_e k = surf_usm_h%k(m) tend(k,j,i) = tend(k,j,i) + g / pt_reference & * surf_usm_h%shf(m) ENDDO ENDIF IF ( use_top_fluxes ) THEN surf_s = surf_def_h(2)%start_index(j,i) surf_e = surf_def_h(2)%end_index(j,i) DO m = surf_s, surf_e k = surf_def_h(2)%k(m) tend(k,j,i) = tend(k,j,i) + g / pt_reference * & surf_def_h(2)%shf(m) ENDDO ENDIF ENDDO ENDIF ENDIF ELSE DO j = nys, nyn DO k = nzb+1, nzt ! !-- Flag 9 is used to mask top fluxes, flag 30 to mask !-- surface fluxes IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN k1 = 1.0_wp + 0.61_wp * q(k,j,i) k2 = 0.61_wp * pt(k,j,i) tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * & g / vpt(k,j,i) * & ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + & k2 * ( q(k+1,j,i) - q(k-1,j,i) ) & ) * dd2zu(k) * & MERGE( 1.0_wp, 0.0_wp, & BTEST( wall_flags_0(k,j,i), 30 ) & ) * & MERGE( 1.0_wp, 0.0_wp, & BTEST( wall_flags_0(k,j,i), 9 ) & ) ELSE IF ( cloud_physics ) THEN IF ( ql(k,j,i) == 0.0_wp ) THEN k1 = 1.0_wp + 0.61_wp * q(k,j,i) k2 = 0.61_wp * pt(k,j,i) ELSE theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i) temp = theta * t_d_pt(k) k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * & ( q(k,j,i) - ql(k,j,i) ) * & ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) / & ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * & ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) ) k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp ) ENDIF tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * & g / vpt(k,j,i) * & ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + & k2 * ( q(k+1,j,i) - q(k-1,j,i) ) & ) * dd2zu(k) * & MERGE( 1.0_wp, 0.0_wp, & BTEST( wall_flags_0(k,j,i), 30 ) & ) * & MERGE( 1.0_wp, 0.0_wp, & BTEST( wall_flags_0(k,j,i), 9 ) & ) ELSE IF ( cloud_droplets ) THEN k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i) k2 = 0.61_wp * pt(k,j,i) tend(k,j,i) = tend(k,j,i) - & kh(k,j,i) * g / vpt(k,j,i) * & ( k1 * ( pt(k+1,j,i)- pt(k-1,j,i) ) + & k2 * ( q(k+1,j,i) - q(k-1,j,i) ) - & pt(k,j,i) * ( ql(k+1,j,i) - & ql(k-1,j,i) ) ) * dd2zu(k) * & MERGE( 1.0_wp, 0.0_wp, & BTEST( wall_flags_0(k,j,i), 30 ) & ) * & MERGE( 1.0_wp, 0.0_wp, & BTEST( wall_flags_0(k,j,i), 9 ) & ) ENDIF ENDDO ENDDO IF ( use_surface_fluxes ) THEN DO j = nys, nyn ! !-- Treat horizontal default surfaces, up- and downward-facing DO l = 0, 1 surf_s = surf_def_h(l)%start_index(j,i) surf_e = surf_def_h(l)%end_index(j,i) DO m = surf_s, surf_e k = surf_def_h(l)%k(m) IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN k1 = 1.0_wp + 0.61_wp * q(k,j,i) k2 = 0.61_wp * pt(k,j,i) ELSE IF ( cloud_physics ) THEN IF ( ql(k,j,i) == 0.0_wp ) THEN k1 = 1.0_wp + 0.61_wp * q(k,j,i) k2 = 0.61_wp * pt(k,j,i) ELSE theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i) temp = theta * t_d_pt(k) k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * & ( q(k,j,i) - ql(k,j,i) ) * & ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) / & ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * & ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) ) k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp ) ENDIF ELSE IF ( cloud_droplets ) THEN k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i) k2 = 0.61_wp * pt(k,j,i) ENDIF tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * & ( k1 * surf_def_h(l)%shf(m) + & k2 * surf_def_h(l)%qsws(m) ) ENDDO ENDDO ! !-- Treat horizontal natural surfaces surf_s = surf_lsm_h%start_index(j,i) surf_e = surf_lsm_h%end_index(j,i) DO m = surf_s, surf_e k = surf_lsm_h%k(m) IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN k1 = 1.0_wp + 0.61_wp * q(k,j,i) k2 = 0.61_wp * pt(k,j,i) ELSE IF ( cloud_physics ) THEN IF ( ql(k,j,i) == 0.0_wp ) THEN k1 = 1.0_wp + 0.61_wp * q(k,j,i) k2 = 0.61_wp * pt(k,j,i) ELSE theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i) temp = theta * t_d_pt(k) k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * & ( q(k,j,i) - ql(k,j,i) ) * & ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) / & ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * & ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) ) k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp ) ENDIF ELSE IF ( cloud_droplets ) THEN k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i) k2 = 0.61_wp * pt(k,j,i) ENDIF tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * & ( k1 * surf_lsm_h%shf(m) + & k2 * surf_lsm_h%qsws(m) ) ENDDO ! !-- Treat horizontal urban surfaces surf_s = surf_usm_h%start_index(j,i) surf_e = surf_usm_h%end_index(j,i) DO m = surf_s, surf_e k = surf_lsm_h%k(m) IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN k1 = 1.0_wp + 0.61_wp * q(k,j,i) k2 = 0.61_wp * pt(k,j,i) ELSE IF ( cloud_physics ) THEN IF ( ql(k,j,i) == 0.0_wp ) THEN k1 = 1.0_wp + 0.61_wp * q(k,j,i) k2 = 0.61_wp * pt(k,j,i) ELSE theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i) temp = theta * t_d_pt(k) k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * & ( q(k,j,i) - ql(k,j,i) ) * & ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) / & ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * & ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) ) k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp ) ENDIF ELSE IF ( cloud_droplets ) THEN k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i) k2 = 0.61_wp * pt(k,j,i) ENDIF tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * & ( k1 * surf_usm_h%shf(m) + & k2 * surf_usm_h%qsws(m) ) ENDDO ENDDO ENDIF IF ( use_top_fluxes ) THEN DO j = nys, nyn surf_s = surf_def_h(2)%start_index(j,i) surf_e = surf_def_h(2)%end_index(j,i) DO m = surf_s, surf_e k = surf_def_h(2)%k(m) IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN k1 = 1.0_wp + 0.61_wp * q(k,j,i) k2 = 0.61_wp * pt(k,j,i) ELSE IF ( cloud_physics ) THEN IF ( ql(k,j,i) == 0.0_wp ) THEN k1 = 1.0_wp + 0.61_wp * q(k,j,i) k2 = 0.61_wp * pt(k,j,i) ELSE theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i) temp = theta * t_d_pt(k) k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * & ( q(k,j,i) - ql(k,j,i) ) * & ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) / & ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * & ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) ) k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp ) ENDIF ELSE IF ( cloud_droplets ) THEN k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i) k2 = 0.61_wp * pt(k,j,i) ENDIF tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * & ( k1 * surf_def_h(2)%shf(m) + & k2 * surf_def_h(2)%qsws(m) ) ENDDO ENDDO ENDIF ENDIF ENDIF ENDDO END SUBROUTINE production_e !------------------------------------------------------------------------------! ! Description: ! ------------ !> Call for grid point i,j !------------------------------------------------------------------------------! SUBROUTINE production_e_ij( i, j ) USE arrays_3d, & ONLY: ddzw, dd2zu, kh, km, prho, pt, q, ql, tend, u, v, vpt, w USE cloud_parameters, & ONLY: l_d_cp, l_d_r, pt_d_t, t_d_pt USE control_parameters, & ONLY: cloud_droplets, cloud_physics, constant_flux_layer, g, & humidity, kappa, neutral, ocean, pt_reference, & rho_reference, use_single_reference_value, & use_surface_fluxes, use_top_fluxes USE grid_variables, & ONLY: ddx, dx, ddy, dy USE indices, & ONLY: nxl, nxr, nys, nyn, nzb, nzb, nzt, wall_flags_0 USE surface_mod, & ONLY : surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, & surf_usm_v IMPLICIT NONE INTEGER(iwp) :: i !< running index x-direction INTEGER(iwp) :: j !< running index y-direction INTEGER(iwp) :: k !< running index z-direction INTEGER(iwp) :: l !< running index for different surface type orientation INTEGER(iwp) :: m !< running index surface elements INTEGER(iwp) :: surf_e !< end index of surface elements at given i-j position INTEGER(iwp) :: surf_s !< start index of surface elements at given i-j position REAL(wp) :: def !< REAL(wp) :: flag !< flag to mask topography REAL(wp) :: k1 !< REAL(wp) :: k2 !< REAL(wp) :: km_neutral !< diffusion coefficient assuming neutral conditions - used to compute shear production at surfaces REAL(wp) :: theta !< REAL(wp) :: temp !< REAL(wp) :: sign_dir !< sign of wall-tke flux, depending on wall orientation REAL(wp) :: usvs !< momentum flux u"v" REAL(wp) :: vsus !< momentum flux v"u" REAL(wp) :: wsus !< momentum flux w"u" REAL(wp) :: wsvs !< momentum flux w"v" REAL(wp), DIMENSION(nzb+1:nzt) :: dudx !< Gradient of u-component in x-direction REAL(wp), DIMENSION(nzb+1:nzt) :: dudy !< Gradient of u-component in y-direction REAL(wp), DIMENSION(nzb+1:nzt) :: dudz !< Gradient of u-component in z-direction REAL(wp), DIMENSION(nzb+1:nzt) :: dvdx !< Gradient of v-component in x-direction REAL(wp), DIMENSION(nzb+1:nzt) :: dvdy !< Gradient of v-component in y-direction REAL(wp), DIMENSION(nzb+1:nzt) :: dvdz !< Gradient of v-component in z-direction REAL(wp), DIMENSION(nzb+1:nzt) :: dwdx !< Gradient of w-component in x-direction REAL(wp), DIMENSION(nzb+1:nzt) :: dwdy !< Gradient of w-component in y-direction REAL(wp), DIMENSION(nzb+1:nzt) :: dwdz !< Gradient of w-component in z-direction IF ( constant_flux_layer ) THEN ! !-- Calculate TKE production by shear. Calculate gradients at all grid !-- points first, gradients at surface-bounded grid points will be !-- overwritten further below. DO k = nzb+1, nzt dudx(k) = ( u(k,j,i+1) - u(k,j,i) ) * ddx dudy(k) = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - & u(k,j-1,i) - u(k,j-1,i+1) ) * ddy dudz(k) = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - & u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k) dvdx(k) = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - & v(k,j,i-1) - v(k,j+1,i-1) ) * ddx dvdy(k) = ( v(k,j+1,i) - v(k,j,i) ) * ddy dvdz(k) = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - & v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k) dwdx(k) = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - & w(k,j,i-1) - w(k-1,j,i-1) ) * ddx dwdy(k) = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - & w(k,j-1,i) - w(k-1,j-1,i) ) * ddy dwdz(k) = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) ENDDO ! !-- Compute gradients at north- and south-facing surfaces. !-- Note, no vertical natural surfaces so far. DO l = 0, 1 ! !-- Default surfaces surf_s = surf_def_v(l)%start_index(j,i) surf_e = surf_def_v(l)%end_index(j,i) DO m = surf_s, surf_e k = surf_def_v(l)%k(m) usvs = surf_def_v(l)%mom_flux_tke(0,m) wsvs = surf_def_v(l)%mom_flux_tke(1,m) km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp & * 0.5_wp * dy ! !-- -1.0 for right-facing wall, 1.0 for left-facing wall sign_dir = MERGE( 1.0_wp, -1.0_wp, & BTEST( wall_flags_0(k,j-1,i), 0 ) ) dudy(k) = sign_dir * usvs / ( km_neutral + 1E-10_wp ) dwdy(k) = sign_dir * wsvs / ( km_neutral + 1E-10_wp ) ENDDO ! !-- Natural surfaces surf_s = surf_lsm_v(l)%start_index(j,i) surf_e = surf_lsm_v(l)%end_index(j,i) DO m = surf_s, surf_e k = surf_lsm_v(l)%k(m) usvs = surf_lsm_v(l)%mom_flux_tke(0,m) wsvs = surf_lsm_v(l)%mom_flux_tke(1,m) km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp & * 0.5_wp * dy ! !-- -1.0 for right-facing wall, 1.0 for left-facing wall sign_dir = MERGE( 1.0_wp, -1.0_wp, & BTEST( wall_flags_0(k,j-1,i), 0 ) ) dudy(k) = sign_dir * usvs / ( km_neutral + 1E-10_wp ) dwdy(k) = sign_dir * wsvs / ( km_neutral + 1E-10_wp ) ENDDO ! !-- Urban surfaces surf_s = surf_usm_v(l)%start_index(j,i) surf_e = surf_usm_v(l)%end_index(j,i) DO m = surf_s, surf_e k = surf_usm_v(l)%k(m) usvs = surf_usm_v(l)%mom_flux_tke(0,m) wsvs = surf_usm_v(l)%mom_flux_tke(1,m) km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp & * 0.5_wp * dy ! !-- -1.0 for right-facing wall, 1.0 for left-facing wall sign_dir = MERGE( 1.0_wp, -1.0_wp, & BTEST( wall_flags_0(k,j-1,i), 0 ) ) dudy(k) = sign_dir * usvs / ( km_neutral + 1E-10_wp ) dwdy(k) = sign_dir * wsvs / ( km_neutral + 1E-10_wp ) ENDDO ENDDO ! !-- Compute gradients at east- and west-facing walls DO l = 2, 3 ! !-- Default surfaces surf_s = surf_def_v(l)%start_index(j,i) surf_e = surf_def_v(l)%end_index(j,i) DO m = surf_s, surf_e k = surf_def_v(l)%k(m) vsus = surf_def_v(l)%mom_flux_tke(0,m) wsus = surf_def_v(l)%mom_flux_tke(1,m) km_neutral = kappa * ( vsus**2 + wsus**2 )**0.25_wp & * 0.5_wp * dx ! !-- -1.0 for right-facing wall, 1.0 for left-facing wall sign_dir = MERGE( 1.0_wp, -1.0_wp, & BTEST( wall_flags_0(k,j,i-1), 0 ) ) dvdx(k) = sign_dir * vsus / ( km_neutral + 1E-10_wp ) dwdx(k) = sign_dir * wsus / ( km_neutral + 1E-10_wp ) ENDDO ! !-- Natural surfaces surf_s = surf_lsm_v(l)%start_index(j,i) surf_e = surf_lsm_v(l)%end_index(j,i) DO m = surf_s, surf_e k = surf_lsm_v(l)%k(m) vsus = surf_lsm_v(l)%mom_flux_tke(0,m) wsus = surf_lsm_v(l)%mom_flux_tke(1,m) km_neutral = kappa * ( vsus**2 + wsus**2 )**0.25_wp & * 0.5_wp * dx ! !-- -1.0 for right-facing wall, 1.0 for left-facing wall sign_dir = MERGE( 1.0_wp, -1.0_wp, & BTEST( wall_flags_0(k,j,i-1), 0 ) ) dvdx(k) = sign_dir * vsus / ( km_neutral + 1E-10_wp ) dwdx(k) = sign_dir * wsus / ( km_neutral + 1E-10_wp ) ENDDO ! !-- Urban surfaces surf_s = surf_usm_v(l)%start_index(j,i) surf_e = surf_usm_v(l)%end_index(j,i) DO m = surf_s, surf_e k = surf_usm_v(l)%k(m) vsus = surf_usm_v(l)%mom_flux_tke(0,m) wsus = surf_usm_v(l)%mom_flux_tke(1,m) km_neutral = kappa * ( vsus**2 + wsus**2 )**0.25_wp & * 0.5_wp * dx ! !-- -1.0 for right-facing wall, 1.0 for left-facing wall sign_dir = MERGE( 1.0_wp, -1.0_wp, & BTEST( wall_flags_0(k,j,i-1), 0 ) ) dvdx(k) = sign_dir * vsus / ( km_neutral + 1E-10_wp ) dwdx(k) = sign_dir * wsus / ( km_neutral + 1E-10_wp ) ENDDO ENDDO ! !-- Compute gradients at upward-facing walls, first for !-- non-natural default surfaces surf_s = surf_def_h(0)%start_index(j,i) surf_e = surf_def_h(0)%end_index(j,i) DO m = surf_s, surf_e k = surf_def_h(0)%k(m) ! !-- Please note, actually, an interpolation of u_0 and v_0 !-- onto the grid center would be required. However, this !-- would require several data transfers between 2D-grid and !-- wall type. The effect of this missing interpolation is !-- negligible. (See also production_e_init). dudz(k) = ( u(k+1,j,i) - surf_def_h(0)%u_0(m) ) * dd2zu(k) dvdz(k) = ( v(k+1,j,i) - surf_def_h(0)%v_0(m) ) * dd2zu(k) ENDDO ! !-- Natural surfaces surf_s = surf_lsm_h%start_index(j,i) surf_e = surf_lsm_h%end_index(j,i) DO m = surf_s, surf_e k = surf_lsm_h%k(m) ! !-- Please note, actually, an interpolation of u_0 and v_0 !-- onto the grid center would be required. However, this !-- would require several data transfers between 2D-grid and !-- wall type. The effect of this missing interpolation is !-- negligible. (See also production_e_init). dudz(k) = ( u(k+1,j,i) - surf_lsm_h%u_0(m) ) * dd2zu(k) dvdz(k) = ( v(k+1,j,i) - surf_lsm_h%v_0(m) ) * dd2zu(k) ENDDO ! !-- Urban surfaces surf_s = surf_usm_h%start_index(j,i) surf_e = surf_usm_h%end_index(j,i) DO m = surf_s, surf_e k = surf_usm_h%k(m) ! !-- Please note, actually, an interpolation of u_0 and v_0 !-- onto the grid center would be required. However, this !-- would require several data transfers between 2D-grid and !-- wall type. The effect of this missing interpolation is !-- negligible. (See also production_e_init). dudz(k) = ( u(k+1,j,i) - surf_usm_h%u_0(m) ) * dd2zu(k) dvdz(k) = ( v(k+1,j,i) - surf_usm_h%v_0(m) ) * dd2zu(k) ENDDO ! !-- Compute gradients at downward-facing walls, only for !-- non-natural default surfaces surf_s = surf_def_h(1)%start_index(j,i) surf_e = surf_def_h(1)%end_index(j,i) DO m = surf_s, surf_e k = surf_def_h(1)%k(m) ! !-- Please note, actually, an interpolation of u_0 and v_0 !-- onto the grid center would be required. However, this !-- would require several data transfers between 2D-grid and !-- wall type. The effect of this missing interpolation is !-- negligible. (See also production_e_init). dudz(k) = ( surf_def_h(1)%u_0(m) - u(k-1,j,i) ) * dd2zu(k) dvdz(k) = ( surf_def_h(1)%v_0(m) - v(k-1,j,i) ) * dd2zu(k) ENDDO DO k = nzb+1, nzt def = 2.0_wp * ( dudx(k)**2 + dvdy(k)**2 + dwdz(k)**2 ) + & dudy(k)**2 + dvdx(k)**2 + dwdx(k)**2 + & dwdy(k)**2 + dudz(k)**2 + dvdz(k)**2 + & 2.0_wp * ( dvdx(k)*dudy(k) + dwdx(k)*dudz(k) + dwdy(k)*dvdz(k) ) IF ( def < 0.0_wp ) def = 0.0_wp flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def * flag ENDDO ELSEIF ( use_surface_fluxes ) THEN ! !-- Calculate TKE production by shear. Here, no additional !-- wall-bounded code is considered. !-- Why? DO k = nzb+1, nzt dudx(k) = ( u(k,j,i+1) - u(k,j,i) ) * ddx dudy(k) = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - & u(k,j-1,i) - u(k,j-1,i+1) ) * ddy dudz(k) = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - & u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k) dvdx(k) = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - & v(k,j,i-1) - v(k,j+1,i-1) ) * ddx dvdy(k) = ( v(k,j+1,i) - v(k,j,i) ) * ddy dvdz(k) = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - & v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k) dwdx(k) = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - & w(k,j,i-1) - w(k-1,j,i-1) ) * ddx dwdy(k) = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - & w(k,j-1,i) - w(k-1,j-1,i) ) * ddy dwdz(k) = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) def = 2.0_wp * ( dudx(k)**2 + dvdy(k)**2 + dwdz(k)**2 ) + & dudy(k)**2 + dvdx(k)**2 + dwdx(k)**2 + & dwdy(k)**2 + dudz(k)**2 + dvdz(k)**2 + & 2.0_wp * ( dvdx(k)*dudy(k) + dwdx(k)*dudz(k) + dwdy(k)*dvdz(k) ) IF ( def < 0.0_wp ) def = 0.0_wp flag = MERGE( 1.0_wp, 0.0_wp, & BTEST( wall_flags_0(k,j,i), 29 ) ) tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def * flag ENDDO ENDIF ! !-- If required, calculate TKE production by buoyancy IF ( .NOT. neutral ) THEN IF ( .NOT. humidity ) THEN IF ( use_single_reference_value ) THEN IF ( ocean ) THEN ! !-- So far in the ocean no special treatment of density flux in !-- the bottom and top surface layer DO k = nzb+1, nzt tend(k,j,i) = tend(k,j,i) + & kh(k,j,i) * g / rho_reference * & ( prho(k+1,j,i) - prho(k-1,j,i) ) * & dd2zu(k) * & MERGE( 1.0_wp, 0.0_wp, & BTEST( wall_flags_0(k,j,i), 0 ) & ) ENDDO ELSE DO k = nzb+1, nzt ! !-- Flag 9 is used to mask top fluxes, flag 30 to mask !-- surface fluxes tend(k,j,i) = tend(k,j,i) - & kh(k,j,i) * g / pt_reference * & ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k) * & MERGE( 1.0_wp, 0.0_wp, & BTEST( wall_flags_0(k,j,i), 30 ) & ) * & MERGE( 1.0_wp, 0.0_wp, & BTEST( wall_flags_0(k,j,i), 9 ) & ) ENDDO IF ( use_surface_fluxes ) THEN ! !-- Default surfaces, up- and downward-facing DO l = 0, 1 surf_s = surf_def_h(l)%start_index(j,i) surf_e = surf_def_h(l)%end_index(j,i) DO m = surf_s, surf_e k = surf_def_h(l)%k(m) tend(k,j,i) = tend(k,j,i) + g / pt_reference * & surf_def_h(l)%shf(m) ENDDO ENDDO ! !-- Natural surfaces surf_s = surf_lsm_h%start_index(j,i) surf_e = surf_lsm_h%end_index(j,i) DO m = surf_s, surf_e k = surf_lsm_h%k(m) tend(k,j,i) = tend(k,j,i) + g / pt_reference * & surf_lsm_h%shf(m) ENDDO ! !-- Urban surfaces surf_s = surf_usm_h%start_index(j,i) surf_e = surf_usm_h%end_index(j,i) DO m = surf_s, surf_e k = surf_usm_h%k(m) tend(k,j,i) = tend(k,j,i) + g / pt_reference * & surf_usm_h%shf(m) ENDDO ENDIF IF ( use_top_fluxes ) THEN surf_s = surf_def_h(2)%start_index(j,i) surf_e = surf_def_h(2)%end_index(j,i) DO m = surf_s, surf_e k = surf_def_h(2)%k(m) tend(k,j,i) = tend(k,j,i) + g / pt_reference * & surf_def_h(2)%shf(m) ENDDO ENDIF ENDIF ELSE IF ( ocean ) THEN ! !-- So far in the ocean no special treatment of density flux in !-- the bottom and top surface layer DO k = nzb+1, nzt tend(k,j,i) = tend(k,j,i) + & kh(k,j,i) * g / prho(k,j,i) * & ( prho(k+1,j,i) - prho(k-1,j,i) ) * & dd2zu(k) * & MERGE( 1.0_wp, 0.0_wp, & BTEST( wall_flags_0(k,j,i), 0 ) & ) ENDDO ELSE DO k = nzb+1, nzt ! !-- Flag 9 is used to mask top fluxes, flag 30 to mask !-- surface fluxes tend(k,j,i) = tend(k,j,i) - & kh(k,j,i) * g / pt(k,j,i) * & ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k) * & MERGE( 1.0_wp, 0.0_wp, & BTEST( wall_flags_0(k,j,i), 30 ) & ) * & MERGE( 1.0_wp, 0.0_wp, & BTEST( wall_flags_0(k,j,i), 9 ) & ) ENDDO IF ( use_surface_fluxes ) THEN ! !-- Default surfaces, up- and downward-facing DO l = 0, 1 surf_s = surf_def_h(l)%start_index(j,i) surf_e = surf_def_h(l)%end_index(j,i) DO m = surf_s, surf_e k = surf_def_h(l)%k(m) tend(k,j,i) = tend(k,j,i) + g / pt_reference & * surf_def_h(l)%shf(m) ENDDO ENDDO ! !-- Natural surfaces surf_s = surf_lsm_h%start_index(j,i) surf_e = surf_lsm_h%end_index(j,i) DO m = surf_s, surf_e k = surf_lsm_h%k(m) tend(k,j,i) = tend(k,j,i) + g / pt_reference & * surf_lsm_h%shf(m) ENDDO ! !-- Urban surfaces surf_s = surf_usm_h%start_index(j,i) surf_e = surf_usm_h%end_index(j,i) DO m = surf_s, surf_e k = surf_usm_h%k(m) tend(k,j,i) = tend(k,j,i) + g / pt_reference & * surf_usm_h%shf(m) ENDDO ENDIF IF ( use_top_fluxes ) THEN surf_s = surf_def_h(2)%start_index(j,i) surf_e = surf_def_h(2)%end_index(j,i) DO m = surf_s, surf_e k = surf_def_h(2)%k(m) tend(k,j,i) = tend(k,j,i) + g / pt_reference * & surf_def_h(2)%shf(m) ENDDO ENDIF ENDIF ENDIF ELSE DO k = nzb+1, nzt ! !-- Flag 9 is used to mask top fluxes, flag 30 to mask surface fluxes IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN k1 = 1.0_wp + 0.61_wp * q(k,j,i) k2 = 0.61_wp * pt(k,j,i) tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) * & ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + & k2 * ( q(k+1,j,i) - q(k-1,j,i) ) & ) * dd2zu(k) * & MERGE( 1.0_wp, 0.0_wp, & BTEST( wall_flags_0(k,j,i), 30 ) & ) * & MERGE( 1.0_wp, 0.0_wp, & BTEST( wall_flags_0(k,j,i), 9 ) & ) ELSE IF ( cloud_physics ) THEN IF ( ql(k,j,i) == 0.0_wp ) THEN k1 = 1.0_wp + 0.61_wp * q(k,j,i) k2 = 0.61_wp * pt(k,j,i) ELSE theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i) temp = theta * t_d_pt(k) k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * & ( q(k,j,i) - ql(k,j,i) ) * & ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) / & ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * & ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) ) k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp ) ENDIF tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) * & ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + & k2 * ( q(k+1,j,i) - q(k-1,j,i) ) & ) * dd2zu(k) * & MERGE( 1.0_wp, 0.0_wp, & BTEST( wall_flags_0(k,j,i), 30 ) & ) * & MERGE( 1.0_wp, 0.0_wp, & BTEST( wall_flags_0(k,j,i), 9 ) & ) ELSE IF ( cloud_droplets ) THEN k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i) k2 = 0.61_wp * pt(k,j,i) tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) * & ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + & k2 * ( q(k+1,j,i) - q(k-1,j,i) ) - & pt(k,j,i) * ( ql(k+1,j,i) - & ql(k-1,j,i) ) ) * dd2zu(k)& * MERGE( 1.0_wp, 0.0_wp, & BTEST( wall_flags_0(k,j,i), 30 ) & ) & * MERGE( 1.0_wp, 0.0_wp, & BTEST( wall_flags_0(k,j,i), 9 ) & ) ENDIF ENDDO IF ( use_surface_fluxes ) THEN ! !-- Treat horizontal default surfaces, up- and downward-facing DO l = 0, 1 surf_s = surf_def_h(l)%start_index(j,i) surf_e = surf_def_h(l)%end_index(j,i) DO m = surf_s, surf_e k = surf_def_h(l)%k(m) IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN k1 = 1.0_wp + 0.61_wp * q(k,j,i) k2 = 0.61_wp * pt(k,j,i) ELSE IF ( cloud_physics ) THEN IF ( ql(k,j,i) == 0.0_wp ) THEN k1 = 1.0_wp + 0.61_wp * q(k,j,i) k2 = 0.61_wp * pt(k,j,i) ELSE theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i) temp = theta * t_d_pt(k) k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * & ( q(k,j,i) - ql(k,j,i) ) * & ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) / & ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * & ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) ) k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp ) ENDIF ELSE IF ( cloud_droplets ) THEN k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i) k2 = 0.61_wp * pt(k,j,i) ENDIF tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * & ( k1 * surf_def_h(l)%shf(m) + & k2 * surf_def_h(l)%qsws(m) ) ENDDO ENDDO ! !-- Treat horizontal natural surfaces surf_s = surf_lsm_h%start_index(j,i) surf_e = surf_lsm_h%end_index(j,i) DO m = surf_s, surf_e k = surf_lsm_h%k(m) IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN k1 = 1.0_wp + 0.61_wp * q(k,j,i) k2 = 0.61_wp * pt(k,j,i) ELSE IF ( cloud_physics ) THEN IF ( ql(k,j,i) == 0.0_wp ) THEN k1 = 1.0_wp + 0.61_wp * q(k,j,i) k2 = 0.61_wp * pt(k,j,i) ELSE theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i) temp = theta * t_d_pt(k) k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * & ( q(k,j,i) - ql(k,j,i) ) * & ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) / & ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * & ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) ) k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp ) ENDIF ELSE IF ( cloud_droplets ) THEN k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i) k2 = 0.61_wp * pt(k,j,i) ENDIF tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * & ( k1 * surf_lsm_h%shf(m) + & k2 * surf_lsm_h%qsws(m) ) ENDDO ! !-- Treat horizontal urban surfaces surf_s = surf_usm_h%start_index(j,i) surf_e = surf_usm_h%end_index(j,i) DO m = surf_s, surf_e k = surf_usm_h%k(m) IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN k1 = 1.0_wp + 0.61_wp * q(k,j,i) k2 = 0.61_wp * pt(k,j,i) ELSE IF ( cloud_physics ) THEN IF ( ql(k,j,i) == 0.0_wp ) THEN k1 = 1.0_wp + 0.61_wp * q(k,j,i) k2 = 0.61_wp * pt(k,j,i) ELSE theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i) temp = theta * t_d_pt(k) k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * & ( q(k,j,i) - ql(k,j,i) ) * & ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) / & ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * & ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) ) k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp ) ENDIF ELSE IF ( cloud_droplets ) THEN k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i) k2 = 0.61_wp * pt(k,j,i) ENDIF tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * & ( k1 * surf_usm_h%shf(m) + & k2 * surf_usm_h%qsws(m) ) ENDDO ENDIF IF ( use_top_fluxes ) THEN surf_s = surf_def_h(2)%start_index(j,i) surf_e = surf_def_h(2)%end_index(j,i) DO m = surf_s, surf_e k = surf_def_h(2)%k(m) IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN k1 = 1.0_wp + 0.61_wp * q(k,j,i) k2 = 0.61_wp * pt(k,j,i) ELSE IF ( cloud_physics ) THEN IF ( ql(k,j,i) == 0.0_wp ) THEN k1 = 1.0_wp + 0.61_wp * q(k,j,i) k2 = 0.61_wp * pt(k,j,i) ELSE theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i) temp = theta * t_d_pt(k) k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * & ( q(k,j,i) - ql(k,j,i) ) * & ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) / & ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * & ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) ) k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp ) ENDIF ELSE IF ( cloud_droplets ) THEN k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i) k2 = 0.61_wp * pt(k,j,i) ENDIF tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * & ( k1* surf_def_h(2)%shf(m) + & k2 * surf_def_h(2)%qsws(m) ) ENDDO ENDIF ENDIF ENDIF END SUBROUTINE production_e_ij !------------------------------------------------------------------------------! ! Description: ! ------------ !> @todo Missing subroutine description. !------------------------------------------------------------------------------! SUBROUTINE production_e_init USE arrays_3d, & ONLY: kh, km, u, v, zu USE control_parameters, & ONLY: constant_flux_layer, kappa USE indices, & ONLY: nxl, nxlg, nxr, nxrg, nys, nysg, nyn, nyng, nzb_u_inner, & nzb_v_inner USE surface_mod, & ONLY : surf_def_h, surf_def_v, surf_lsm_h, surf_usm_h IMPLICIT NONE INTEGER(iwp) :: i !< grid index x-direction INTEGER(iwp) :: j !< grid index y-direction INTEGER(iwp) :: k !< grid index z-direction INTEGER(iwp) :: l !< running index surface type (up- or downward-facing) INTEGER(iwp) :: m !< running index surface elements IF ( constant_flux_layer ) THEN ! !-- Calculate a virtual velocity at the surface in a way that the !-- vertical velocity gradient at k = 1 (u(k+1)-u_0) matches the !-- Prandtl law (-w'u'/km). This gradient is used in the TKE shear !-- production term at k=1 (see production_e_ij). !-- The velocity gradient has to be limited in case of too small km !-- (otherwise the timestep may be significantly reduced by large !-- surface winds). !-- not available in case of non-cyclic boundary conditions. !-- WARNING: the exact analytical solution would require the determination !-- of the eddy diffusivity by km = u* * kappa * zp / phi_m. !-- Default surfaces, upward-facing !$OMP PARALLEL DO PRIVATE(i,j,k,m) DO m = 1, surf_def_h(0)%ns i = surf_def_h(0)%i(m) j = surf_def_h(0)%j(m) k = surf_def_h(0)%k(m) ! !-- Note, calculatione of u_0 and v_0 is not fully accurate, as u/v !-- and km are not on the same grid. Actually, a further !-- interpolation of km onto the u/v-grid is necessary. However, the !-- effect of this error is negligible. surf_def_h(0)%u_0(m) = u(k+1,j,i) + surf_def_h(0)%usws(m) * & ( zu(k+1) - zu(k-1) ) / & ( km(k,j,i) + 1.0E-20_wp ) surf_def_h(0)%v_0(m) = v(k+1,j,i) + surf_def_h(0)%vsws(m) * & ( zu(k+1) - zu(k-1) ) / & ( km(k,j,i) + 1.0E-20_wp ) IF ( ABS( u(k+1,j,i) - surf_def_h(0)%u_0(m) ) > & ABS( u(k+1,j,i) - u(k-1,j,i) ) & ) surf_def_h(0)%u_0(m) = u(k-1,j,i) IF ( ABS( v(k+1,j,i) - surf_def_h(0)%v_0(m) ) > & ABS( v(k+1,j,i) - v(k-1,j,i) ) & ) surf_def_h(0)%v_0(m) = v(k-1,j,i) ENDDO ! !-- Default surfaces, downward-facing !$OMP PARALLEL DO PRIVATE(i,j,k,m) DO m = 1, surf_def_h(1)%ns i = surf_def_h(1)%i(m) j = surf_def_h(1)%j(m) k = surf_def_h(1)%k(m) ! !-- Note, calculatione of u_0 and v_0 is not fully accurate, as u/v !-- and km are not on the same grid. Actually, a further !-- interpolation of km onto the u/v-grid is necessary. However, the !-- effect of this error is negligible. !-- In case of downward-facing surfaces, gradient is calculated !-- between u_0 and u(k-1). surf_def_h(1)%u_0(m) = u(k-1,j,i) - surf_def_h(1)%usws(m) * & ( zu(k+1) - zu(k-1) ) / & ( km(k,j,i) + 1.0E-20_wp ) surf_def_h(1)%v_0(m) = v(k-1,j,i) - surf_def_h(1)%vsws(m) * & ( zu(k+1) - zu(k-1) ) / & ( km(k,j,i) + 1.0E-20_wp ) IF ( ABS( surf_def_h(1)%u_0(m) - u(k-1,j,i) ) > & ABS( u(k+1,j,i) - u(k-1,j,i) ) & ) surf_def_h(1)%u_0(m) = u(k+1,j,i) IF ( ABS( surf_def_h(1)%v_0(m) - v(k-1,j,i) ) > & ABS( v(k+1,j,i) - v(k-1,j,i) ) & ) surf_def_h(1)%v_0(m) = v(k+1,j,i) ENDDO ! !-- Natural surfaces, upward-facing !$OMP PARALLEL DO PRIVATE(i,j,k,m) DO m = 1, surf_lsm_h%ns i = surf_lsm_h%i(m) j = surf_lsm_h%j(m) k = surf_lsm_h%k(m) ! !-- Note, calculatione of u_0 and v_0 is not fully accurate, as u/v !-- and km are not on the same grid. Actually, a further !-- interpolation of km onto the u/v-grid is necessary. However, the !-- effect of this error is negligible. surf_lsm_h%u_0(m) = u(k+1,j,i) + surf_lsm_h%usws(m) * & ( zu(k+1) - zu(k-1) ) / & ( km(k,j,i) + 1.0E-20_wp ) surf_lsm_h%v_0(m) = v(k+1,j,i) + surf_lsm_h%vsws(m) * & ( zu(k+1) - zu(k-1) ) / & ( km(k,j,i) + 1.0E-20_wp ) IF ( ABS( u(k+1,j,i) - surf_lsm_h%u_0(m) ) > & ABS( u(k+1,j,i) - u(k-1,j,i) ) & ) surf_lsm_h%u_0(m) = u(k-1,j,i) IF ( ABS( v(k+1,j,i) - surf_lsm_h%v_0(m) ) > & ABS( v(k+1,j,i) - v(k-1,j,i) ) & ) surf_lsm_h%v_0(m) = v(k-1,j,i) ENDDO ! !-- Urban surfaces, upward-facing !$OMP PARALLEL DO PRIVATE(i,j,k,m) DO m = 1, surf_usm_h%ns i = surf_usm_h%i(m) j = surf_usm_h%j(m) k = surf_usm_h%k(m) ! !-- Note, calculatione of u_0 and v_0 is not fully accurate, as u/v !-- and km are not on the same grid. Actually, a further !-- interpolation of km onto the u/v-grid is necessary. However, the !-- effect of this error is negligible. surf_usm_h%u_0(m) = u(k+1,j,i) + surf_usm_h%usws(m) * & ( zu(k+1) - zu(k-1) ) / & ( km(k,j,i) + 1.0E-20_wp ) surf_usm_h%v_0(m) = v(k+1,j,i) + surf_usm_h%vsws(m) * & ( zu(k+1) - zu(k-1) ) / & ( km(k,j,i) + 1.0E-20_wp ) IF ( ABS( u(k+1,j,i) - surf_usm_h%u_0(m) ) > & ABS( u(k+1,j,i) - u(k-1,j,i) ) & ) surf_usm_h%u_0(m) = u(k-1,j,i) IF ( ABS( v(k+1,j,i) - surf_usm_h%v_0(m) ) > & ABS( v(k+1,j,i) - v(k-1,j,i) ) & ) surf_usm_h%v_0(m) = v(k-1,j,i) ENDDO ENDIF END SUBROUTINE production_e_init END MODULE production_e_mod