Home

Context Navigation

source: palm/trunk/SOURCE/production_e.f90 @ 2481

Last change on this file since 2481 was 2478, checked in by suehring, 7 years ago
Bugfixes concerning top fluxes and TKE production
Property svn:keywords set to `Id`
File size: 85.3 KB

Rev	Line
[1873]	1	!> @file production_e.f90
[2000]	2	!------------------------------------------------------------------------------!
[1036]	3	! This file is part of PALM.
	4	!
[2000]	5	! PALM is free software: you can redistribute it and/or modify it under the
	6	! terms of the GNU General Public License as published by the Free Software
	7	! Foundation, either version 3 of the License, or (at your option) any later
	8	! version.
[1036]	9	!
	10	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
	11	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
	12	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
	13	!
	14	! You should have received a copy of the GNU General Public License along with
	15	! PALM. If not, see <http://www.gnu.org/licenses/>.
	16	!
[2101]	17	! Copyright 1997-2017 Leibniz Universitaet Hannover
[2000]	18	!------------------------------------------------------------------------------!
[1036]	19	!
[484]	20	! Current revisions:
[1]	21	! -----------------
[2233]	22	!
	23	!
[1321]	24	! Former revisions:
	25	! -----------------
	26	! $Id: production_e.f90 2478 2017-09-18 13:37:24Z raasch $
[2478]	27	! Bugfix, consider case where no constant-flux layer and no surfaces fluxes
	28	! are used
	29	!
	30	! 2329 2017-08-03 14:24:56Z knoop
[2329]	31	! Bugfix: added division by density as kinematic fluxes are needed
	32	!
	33	! 2233 2017-05-30 18:08:54Z suehring
[1321]	34	!
[2233]	35	! 2232 2017-05-30 17:47:52Z suehring
	36	! Adjustments to new surface concept
	37	!
[2127]	38	! 2126 2017-01-20 15:54:21Z raasch
	39	! density in ocean case replaced by potential density
	40	!
[2119]	41	! 2118 2017-01-17 16:38:49Z raasch
	42	! OpenACC version of subroutine removed
	43	!
[2032]	44	! 2031 2016-10-21 15:11:58Z knoop
	45	! renamed variable rho to rho_ocean
	46	!
[2001]	47	! 2000 2016-08-20 18:09:15Z knoop
	48	! Forced header and separation lines into 80 columns
	49	!
[1874]	50	! 1873 2016-04-18 14:50:06Z maronga
	51	! Module renamed (removed _mod)
	52	!
	53	!
[1851]	54	! 1850 2016-04-08 13:29:27Z maronga
	55	! Module renamed
	56	!
	57	!
[1692]	58	! 1691 2015-10-26 16:17:44Z maronga
	59	! Renamed prandtl_layer to constant_flux_layer.
	60	!
[1683]	61	! 1682 2015-10-07 23:56:08Z knoop
	62	! Code annotations made doxygen readable
	63	!
[1375]	64	! 1374 2014-04-25 12:55:07Z raasch
	65	! nzb_s_outer removed from acc-present-list
	66	!
[1354]	67	! 1353 2014-04-08 15:21:23Z heinze
	68	! REAL constants provided with KIND-attribute
	69	!
[1343]	70	! 1342 2014-03-26 17:04:47Z kanani
	71	! REAL constants defined as wp-kind
	72	!
[1321]	73	! 1320 2014-03-20 08:40:49Z raasch
[1320]	74	! ONLY-attribute added to USE-statements,
	75	! kind-parameters added to all INTEGER and REAL declaration statements,
	76	! kinds are defined in new module kinds,
	77	! old module precision_kind is removed,
	78	! revision history before 2012 removed,
	79	! comment fields (!:) to be used for variable explanations added to
	80	! all variable declaration statements
[110]	81	!
[1258]	82	! 1257 2013-11-08 15:18:40Z raasch
	83	! openacc loop and loop vector clauses removed, declare create moved after
	84	! the FORTRAN declaration statement
	85	!
[1182]	86	! 1179 2013-06-14 05:57:58Z raasch
	87	! use_reference renamed use_single_reference_value
	88	!
[1132]	89	! 1128 2013-04-12 06:19:32Z raasch
	90	! loop index bounds in accelerator version replaced by i_left, i_right, j_south,
	91	! j_north
	92	!
[1037]	93	! 1036 2012-10-22 13:43:42Z raasch
	94	! code put under GPL (PALM 3.9)
	95	!
[1017]	96	! 1015 2012-09-27 09:23:24Z raasch
	97	! accelerator version (*_acc) added
	98	!
[1008]	99	! 1007 2012-09-19 14:30:36Z franke
	100	! Bugfix: calculation of buoyancy production has to consider the liquid water
	101	! mixing ratio in case of cloud droplets
	102	!
[941]	103	! 940 2012-07-09 14:31:00Z raasch
	104	! TKE production by buoyancy can be switched off in case of runs with pure
	105	! neutral stratification
	106	!
[1]	107	! Revision 1.1 1997/09/19 07:45:35 raasch
	108	! Initial revision
	109	!
	110	!
	111	! Description:
	112	! ------------
[1682]	113	!> Production terms (shear + buoyancy) of the TKE.
[1691]	114	!> @warning The case with constant_flux_layer = F and use_surface_fluxes = T is
[1682]	115	!> not considered well!
[1]	116	!------------------------------------------------------------------------------!
[1682]	117	MODULE production_e_mod
[1]	118
[1320]	119	USE kinds
	120
[1]	121	PRIVATE
[2118]	122	PUBLIC production_e, production_e_init
[56]	123
[1]	124	INTERFACE production_e
	125	MODULE PROCEDURE production_e
	126	MODULE PROCEDURE production_e_ij
	127	END INTERFACE production_e
	128
	129	INTERFACE production_e_init
	130	MODULE PROCEDURE production_e_init
	131	END INTERFACE production_e_init
	132
	133	CONTAINS
	134
	135
	136	!------------------------------------------------------------------------------!
[1682]	137	! Description:
	138	! ------------
	139	!> Call for all grid points
[1]	140	!------------------------------------------------------------------------------!
	141	SUBROUTINE production_e
	142
[1320]	143	USE arrays_3d, &
[2329]	144	ONLY: ddzw, dd2zu, drho_air_zw, kh, km, prho, pt, q, ql, tend, u, &
	145	v, vpt, w
[1]	146
[1320]	147	USE cloud_parameters, &
	148	ONLY: l_d_cp, l_d_r, pt_d_t, t_d_pt
	149
	150	USE control_parameters, &
[1691]	151	ONLY: cloud_droplets, cloud_physics, constant_flux_layer, g, &
	152	humidity, kappa, neutral, ocean, pt_reference, &
	153	rho_reference, use_single_reference_value, &
	154	use_surface_fluxes, use_top_fluxes
[1320]	155
	156	USE grid_variables, &
[2232]	157	ONLY: ddx, dx, ddy, dy
[1320]	158
	159	USE indices, &
[2232]	160	ONLY: nxl, nxr, nys, nyn, nzb, nzt, wall_flags_0
[1320]	161
[2232]	162	USE surface_mod, &
	163	ONLY : surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, &
	164	surf_usm_v
	165
[1]	166	IMPLICIT NONE
	167
[2232]	168	INTEGER(iwp) :: i !< running index x-direction
	169	INTEGER(iwp) :: j !< running index y-direction
	170	INTEGER(iwp) :: k !< running index z-direction
	171	INTEGER(iwp) :: l !< running index for different surface type orientation
	172	INTEGER(iwp) :: m !< running index surface elements
	173	INTEGER(iwp) :: surf_e !< end index of surface elements at given i-j position
	174	INTEGER(iwp) :: surf_s !< start index of surface elements at given i-j position
[1]	175
[2232]	176	REAL(wp) :: def !<
	177	REAL(wp) :: flag !< flag to mask topography
[1682]	178	REAL(wp) :: k1 !<
	179	REAL(wp) :: k2 !<
[2232]	180	REAL(wp) :: km_neutral !< diffusion coefficient assuming neutral conditions - used to compute shear production at surfaces
[1682]	181	REAL(wp) :: theta !<
	182	REAL(wp) :: temp !<
[2232]	183	REAL(wp) :: sign_dir !< sign of wall-tke flux, depending on wall orientation
	184	REAL(wp) :: usvs !< momentum flux u"v"
	185	REAL(wp) :: vsus !< momentum flux v"u"
	186	REAL(wp) :: wsus !< momentum flux w"u"
	187	REAL(wp) :: wsvs !< momentum flux w"v"
[1]	188
[2232]	189	REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: dudx !< Gradient of u-component in x-direction
	190	REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: dudy !< Gradient of u-component in y-direction
	191	REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: dudz !< Gradient of u-component in z-direction
	192	REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: dvdx !< Gradient of v-component in x-direction
	193	REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: dvdy !< Gradient of v-component in y-direction
	194	REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: dvdz !< Gradient of v-component in z-direction
	195	REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: dwdx !< Gradient of w-component in x-direction
	196	REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: dwdy !< Gradient of w-component in y-direction
	197	REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: dwdz !< Gradient of w-component in z-direction
[1]	198
[2232]	199	DO i = nxl, nxr
[53]	200
[2232]	201	IF ( constant_flux_layer ) THEN
[940]	202
	203	!
[2232]	204	!-- Calculate TKE production by shear. Calculate gradients at all grid
	205	!-- points first, gradients at surface-bounded grid points will be
	206	!-- overwritten further below.
	207	DO j = nys, nyn
	208	DO k = nzb+1, nzt
[1]	209
[2232]	210	dudx(k,j) = ( u(k,j,i+1) - u(k,j,i) ) * ddx
	211	dudy(k,j) = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
	212	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
	213	dudz(k,j) = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - &
	214	u(k-1,j,i) - u(k-1,j,i+1) ) * &
	215	dd2zu(k)
	216
	217	dvdx(k,j) = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
	218	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
	219	dvdy(k,j) = ( v(k,j+1,i) - v(k,j,i) ) * ddy
	220	dvdz(k,j) = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - &
	221	v(k-1,j,i) - v(k-1,j+1,i) ) * &
	222	dd2zu(k)
[1]	223
[2232]	224	dwdx(k,j) = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
	225	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
	226	dwdy(k,j) = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
	227	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
	228	dwdz(k,j) = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
[1]	229
[2232]	230	ENDDO
[1]	231	ENDDO
	232
	233	!
[55]	234	!-- Position beneath wall
	235	!-- (2) - Will allways be executed.
	236	!-- 'bottom and wall: use u_0,v_0 and wall functions'
[1]	237	DO j = nys, nyn
[1007]	238	!
[2232]	239	!-- Compute gradients at north- and south-facing surfaces.
	240	!-- First, for default surfaces, then for urban surfaces.
	241	!-- Note, so far no natural vertical surfaces implemented
	242	DO l = 0, 1
	243	surf_s = surf_def_v(l)%start_index(j,i)
	244	surf_e = surf_def_v(l)%end_index(j,i)
	245	DO m = surf_s, surf_e
	246	k = surf_def_v(l)%k(m)
	247	usvs = surf_def_v(l)%mom_flux_tke(0,m)
	248	wsvs = surf_def_v(l)%mom_flux_tke(1,m)
	249
	250	km_neutral = kappa * ( usvs2 + wsvs2 )**0.25_wp &
	251	* 0.5_wp * dy
[1007]	252	!
[2232]	253	!-- -1.0 for right-facing wall, 1.0 for left-facing wall
	254	sign_dir = MERGE( 1.0_wp, -1.0_wp, &
	255	BTEST( wall_flags_0(k,j-1,i), 0 ) )
	256	dudy(k,j) = sign_dir * usvs / ( km_neutral + 1E-10_wp )
	257	dwdy(k,j) = sign_dir * wsvs / ( km_neutral + 1E-10_wp )
	258	ENDDO
[1]	259	!
[2232]	260	!-- Natural surfaces
	261	surf_s = surf_lsm_v(l)%start_index(j,i)
	262	surf_e = surf_lsm_v(l)%end_index(j,i)
	263	DO m = surf_s, surf_e
	264	k = surf_lsm_v(l)%k(m)
	265	usvs = surf_lsm_v(l)%mom_flux_tke(0,m)
	266	wsvs = surf_lsm_v(l)%mom_flux_tke(1,m)
	267
	268	km_neutral = kappa * ( usvs2 + wsvs2 )**0.25_wp &
	269	* 0.5_wp * dy
	270	!
	271	!-- -1.0 for right-facing wall, 1.0 for left-facing wall
	272	sign_dir = MERGE( 1.0_wp, -1.0_wp, &
	273	BTEST( wall_flags_0(k,j-1,i), 0 ) )
	274	dudy(k,j) = sign_dir * usvs / ( km_neutral + 1E-10_wp )
	275	dwdy(k,j) = sign_dir * wsvs / ( km_neutral + 1E-10_wp )
	276	ENDDO
	277	!
	278	!-- Urban surfaces
	279	surf_s = surf_usm_v(l)%start_index(j,i)
	280	surf_e = surf_usm_v(l)%end_index(j,i)
	281	DO m = surf_s, surf_e
	282	k = surf_usm_v(l)%k(m)
	283	usvs = surf_usm_v(l)%mom_flux_tke(0,m)
	284	wsvs = surf_usm_v(l)%mom_flux_tke(1,m)
	285
	286	km_neutral = kappa * ( usvs2 + wsvs2 )**0.25_wp &
	287	* 0.5_wp * dy
	288	!
	289	!-- -1.0 for right-facing wall, 1.0 for left-facing wall
	290	sign_dir = MERGE( 1.0_wp, -1.0_wp, &
	291	BTEST( wall_flags_0(k,j-1,i), 0 ) )
	292	dudy(k,j) = sign_dir * usvs / ( km_neutral + 1E-10_wp )
	293	dwdy(k,j) = sign_dir * wsvs / ( km_neutral + 1E-10_wp )
	294	ENDDO
	295	ENDDO
	296	!
	297	!-- Compute gradients at east- and west-facing walls
	298	DO l = 2, 3
	299	surf_s = surf_def_v(l)%start_index(j,i)
	300	surf_e = surf_def_v(l)%end_index(j,i)
	301	DO m = surf_s, surf_e
	302	k = surf_def_v(l)%k(m)
	303	vsus = surf_def_v(l)%mom_flux_tke(0,m)
	304	wsus = surf_def_v(l)%mom_flux_tke(1,m)
[1]	305
[2232]	306	km_neutral = kappa * ( vsus2 + wsus2 )**0.25_wp &
	307	* 0.5_wp * dx
[1007]	308	!
[2232]	309	!-- -1.0 for right-facing wall, 1.0 for left-facing wall
	310	sign_dir = MERGE( 1.0_wp, -1.0_wp, &
	311	BTEST( wall_flags_0(k,j,i-1), 0 ) )
	312	dvdx(k,j) = sign_dir * vsus / ( km_neutral + 1E-10_wp )
	313	dwdx(k,j) = sign_dir * wsus / ( km_neutral + 1E-10_wp )
	314	ENDDO
	315	!
	316	!-- Natural surfaces
	317	surf_s = surf_lsm_v(l)%start_index(j,i)
	318	surf_e = surf_lsm_v(l)%end_index(j,i)
	319	DO m = surf_s, surf_e
	320	k = surf_lsm_v(l)%k(m)
	321	vsus = surf_lsm_v(l)%mom_flux_tke(0,m)
	322	wsus = surf_lsm_v(l)%mom_flux_tke(1,m)
[1]	323
[2232]	324	km_neutral = kappa * ( vsus2 + wsus2 )**0.25_wp &
	325	* 0.5_wp * dx
[1007]	326	!
[2232]	327	!-- -1.0 for right-facing wall, 1.0 for left-facing wall
	328	sign_dir = MERGE( 1.0_wp, -1.0_wp, &
	329	BTEST( wall_flags_0(k,j,i-1), 0 ) )
	330	dvdx(k,j) = sign_dir * vsus / ( km_neutral + 1E-10_wp )
	331	dwdx(k,j) = sign_dir * wsus / ( km_neutral + 1E-10_wp )
	332	ENDDO
	333	!
	334	!-- Urban surfaces
	335	surf_s = surf_usm_v(l)%start_index(j,i)
	336	surf_e = surf_usm_v(l)%end_index(j,i)
	337	DO m = surf_s, surf_e
	338	k = surf_usm_v(l)%k(m)
	339	vsus = surf_usm_v(l)%mom_flux_tke(0,m)
	340	wsus = surf_usm_v(l)%mom_flux_tke(1,m)
[1]	341
[2232]	342	km_neutral = kappa * ( vsus2 + wsus2 )**0.25_wp &
	343	* 0.5_wp * dx
[1]	344	!
[2232]	345	!-- -1.0 for right-facing wall, 1.0 for left-facing wall
	346	sign_dir = MERGE( 1.0_wp, -1.0_wp, &
	347	BTEST( wall_flags_0(k,j,i-1), 0 ) )
	348	dvdx(k,j) = sign_dir * vsus / ( km_neutral + 1E-10_wp )
	349	dwdx(k,j) = sign_dir * wsus / ( km_neutral + 1E-10_wp )
	350	ENDDO
	351	ENDDO
	352	!
	353	!-- Compute gradients at upward-facing surfaces
	354	surf_s = surf_def_h(0)%start_index(j,i)
	355	surf_e = surf_def_h(0)%end_index(j,i)
	356	DO m = surf_s, surf_e
	357	k = surf_def_h(0)%k(m)
	358	!
	359	!-- Please note, actually, an interpolation of u_0 and v_0
	360	!-- onto the grid center would be required. However, this
	361	!-- would require several data transfers between 2D-grid and
	362	!-- wall type. The effect of this missing interpolation is
	363	!-- negligible. (See also production_e_init).
	364	dudz(k,j) = ( u(k+1,j,i) - surf_def_h(0)%u_0(m) ) * dd2zu(k)
	365	dvdz(k,j) = ( v(k+1,j,i) - surf_def_h(0)%v_0(m) ) * dd2zu(k)
	366
	367	ENDDO
	368	!
	369	!-- Natural surfaces
	370	surf_s = surf_lsm_h%start_index(j,i)
	371	surf_e = surf_lsm_h%end_index(j,i)
	372	DO m = surf_s, surf_e
	373	k = surf_lsm_h%k(m)
	374	!
	375	!-- Please note, actually, an interpolation of u_0 and v_0
	376	!-- onto the grid center would be required. However, this
	377	!-- would require several data transfers between 2D-grid and
	378	!-- wall type. The effect of this missing interpolation is
	379	!-- negligible. (See also production_e_init).
	380	dudz(k,j) = ( u(k+1,j,i) - surf_lsm_h%u_0(m) ) * dd2zu(k)
	381	dvdz(k,j) = ( v(k+1,j,i) - surf_lsm_h%v_0(m) ) * dd2zu(k)
	382
	383	ENDDO
	384	!
	385	!-- Urban surfaces
	386	surf_s = surf_usm_h%start_index(j,i)
	387	surf_e = surf_usm_h%end_index(j,i)
	388	DO m = surf_s, surf_e
	389	k = surf_usm_h%k(m)
	390	!
	391	!-- Please note, actually, an interpolation of u_0 and v_0
	392	!-- onto the grid center would be required. However, this
	393	!-- would require several data transfers between 2D-grid and
	394	!-- wall type. The effect of this missing interpolation is
	395	!-- negligible. (See also production_e_init).
	396	dudz(k,j) = ( u(k+1,j,i) - surf_usm_h%u_0(m) ) * dd2zu(k)
	397	dvdz(k,j) = ( v(k+1,j,i) - surf_usm_h%v_0(m) ) * dd2zu(k)
	398
	399	ENDDO
	400	!
	401	!-- Compute gradients at downward-facing walls, only for
	402	!-- non-natural default surfaces
	403	surf_s = surf_def_h(1)%start_index(j,i)
	404	surf_e = surf_def_h(1)%end_index(j,i)
	405	DO m = surf_s, surf_e
	406	k = surf_def_h(1)%k(m)
	407	!
	408	!-- Please note, actually, an interpolation of u_0 and v_0
	409	!-- onto the grid center would be required. However, this
	410	!-- would require several data transfers between 2D-grid and
	411	!-- wall type. The effect of this missing interpolation is
	412	!-- negligible. (See also production_e_init).
	413	dudz(k,j) = ( surf_def_h(1)%u_0(m) - u(k-1,j,i) ) * dd2zu(k)
	414	dvdz(k,j) = ( surf_def_h(1)%v_0(m) - v(k-1,j,i) ) * dd2zu(k)
[1]	415
[2232]	416	ENDDO
[1]	417
	418	ENDDO
	419
	420	DO j = nys, nyn
[2232]	421	DO k = nzb+1, nzt
[1]	422
[2232]	423	def = 2.0_wp * ( dudx(k,j)2 + dvdy(k,j)2 + dwdz(k,j)**2 ) + &
	424	dudy(k,j)2 + dvdx(k,j)2 + dwdx(k,j)**2 + &
	425	dwdy(k,j)2 + dudz(k,j)2 + dvdz(k,j)**2 + &
	426	2.0_wp * ( dvdx(k,j)dudy(k,j) + dwdx(k,j)dudz(k,j) + &
	427	dwdy(k,j)*dvdz(k,j) )
[1]	428
[1342]	429	IF ( def < 0.0_wp ) def = 0.0_wp
[1]	430
[2232]	431	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
[1007]	432
[2232]	433	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def * flag
[1]	434
[2232]	435	ENDDO
[1]	436	ENDDO
	437
[2478]	438	ELSE
[37]	439
	440	DO j = nys, nyn
[2232]	441	!
	442	!-- Calculate TKE production by shear. Here, no additional
	443	!-- wall-bounded code is considered.
	444	!-- Why?
	445	DO k = nzb+1, nzt
[37]	446
[2232]	447	dudx(k,j) = ( u(k,j,i+1) - u(k,j,i) ) * ddx
	448	dudy(k,j) = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
	449	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
	450	dudz(k,j) = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - &
	451	u(k-1,j,i) - u(k-1,j,i+1) ) * &
	452	dd2zu(k)
[37]	453
[2232]	454	dvdx(k,j) = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
	455	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
	456	dvdy(k,j) = ( v(k,j+1,i) - v(k,j,i) ) * ddy
	457	dvdz(k,j) = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - &
	458	v(k-1,j,i) - v(k-1,j+1,i) ) * &
	459	dd2zu(k)
[37]	460
[2232]	461	dwdx(k,j) = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
	462	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
	463	dwdy(k,j) = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
	464	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
	465	dwdz(k,j) = ( w(k,j,i) - w(k-1,j,i) ) * &
	466	ddzw(k)
	467
	468	def = 2.0_wp * ( &
	469	dudx(k,j)2 + dvdy(k,j)2 + dwdz(k,j)**2 &
	470	) + &
	471	dudy(k,j)2 + dvdx(k,j)2 + dwdx(k,j)**2 + &
	472	dwdy(k,j)2 + dudz(k,j)2 + dvdz(k,j)**2 + &
	473	2.0_wp * ( &
	474	dvdx(k,j)dudy(k,j) + dwdx(k,j)dudz(k,j) + &
	475	dwdy(k,j)*dvdz(k,j) &
	476	)
[37]	477
[2232]	478	IF ( def < 0.0_wp ) def = 0.0_wp
[37]	479
[2232]	480	flag = MERGE( 1.0_wp, 0.0_wp, &
	481	BTEST( wall_flags_0(k,j,i), 29 ) )
	482	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def * flag
	483
	484	ENDDO
[37]	485	ENDDO
	486
[1]	487	ENDIF
	488
	489	!
[940]	490	!-- If required, calculate TKE production by buoyancy
	491	IF ( .NOT. neutral ) THEN
[1]	492
[940]	493	IF ( .NOT. humidity ) THEN
[1]	494
[1179]	495	IF ( use_single_reference_value ) THEN
[940]	496
	497	IF ( ocean ) THEN
[97]	498	!
[940]	499	!-- So far in the ocean no special treatment of density flux
	500	!-- in the bottom and top surface layer
	501	DO j = nys, nyn
[2232]	502	DO k = nzb+1, nzt
[2126]	503	tend(k,j,i) = tend(k,j,i) + &
	504	kh(k,j,i) * g / rho_reference * &
	505	( prho(k+1,j,i) - prho(k-1,j,i) ) * &
[2232]	506	dd2zu(k) * &
	507	MERGE( 1.0_wp, 0.0_wp, &
	508	BTEST( wall_flags_0(k,j,i), 0 ) &
	509	)
[940]	510	ENDDO
[97]	511	ENDDO
	512
[940]	513	ELSE
[97]	514
[940]	515	DO j = nys, nyn
[2232]	516	DO k = nzb+1, nzt
	517	!
	518	!-- Flag 9 is used to mask top fluxes, flag 30 to mask
	519	!-- surface fluxes
	520	tend(k,j,i) = tend(k,j,i) - &
	521	kh(k,j,i) * g / pt_reference * &
	522	( pt(k+1,j,i) - pt(k-1,j,i) ) * &
	523	dd2zu(k) * &
	524	MERGE( 1.0_wp, 0.0_wp, &
	525	BTEST( wall_flags_0(k,j,i), 30 ) &
	526	) * &
	527	MERGE( 1.0_wp, 0.0_wp, &
	528	BTEST( wall_flags_0(k,j,i), 9 ) &
	529	)
[940]	530	ENDDO
[97]	531
[940]	532	IF ( use_surface_fluxes ) THEN
[2232]	533	!
	534	!-- Default surfaces, up- and downward-facing
	535	DO l = 0, 1
	536	surf_s = surf_def_h(l)%start_index(j,i)
	537	surf_e = surf_def_h(l)%end_index(j,i)
	538	DO m = surf_s, surf_e
	539	k = surf_def_h(l)%k(m)
	540	tend(k,j,i) = tend(k,j,i) + g / pt_reference &
[2329]	541	* drho_air_zw(k-1) &
[2232]	542	* surf_def_h(l)%shf(m)
	543	ENDDO
	544	ENDDO
	545	!
	546	!-- Natural surfaces
	547	surf_s = surf_lsm_h%start_index(j,i)
	548	surf_e = surf_lsm_h%end_index(j,i)
	549	DO m = surf_s, surf_e
	550	k = surf_lsm_h%k(m)
	551	tend(k,j,i) = tend(k,j,i) + g / pt_reference &
[2329]	552	* drho_air_zw(k-1) &
[2232]	553	* surf_lsm_h%shf(m)
	554	ENDDO
	555	!
	556	!-- Urban surfaces
	557	surf_s = surf_usm_h%start_index(j,i)
	558	surf_e = surf_usm_h%end_index(j,i)
	559	DO m = surf_s, surf_e
	560	k = surf_usm_h%k(m)
	561	tend(k,j,i) = tend(k,j,i) + g / pt_reference &
[2329]	562	* drho_air_zw(k-1) &
[2232]	563	* surf_usm_h%shf(m)
	564	ENDDO
[940]	565	ENDIF
[97]	566
[940]	567	IF ( use_top_fluxes ) THEN
[2232]	568	surf_s = surf_def_h(2)%start_index(j,i)
	569	surf_e = surf_def_h(2)%end_index(j,i)
	570	DO m = surf_s, surf_e
	571	k = surf_def_h(2)%k(m)
[2329]	572	tend(k,j,i) = tend(k,j,i) + g / pt_reference &
	573	* drho_air_zw(k-1) &
	574	* surf_def_h(2)%shf(m)
[2232]	575	ENDDO
[940]	576	ENDIF
	577	ENDDO
[57]	578
[940]	579	ENDIF
[57]	580
[940]	581	ELSE
[1]	582
[940]	583	IF ( ocean ) THEN
[97]	584	!
[940]	585	!-- So far in the ocean no special treatment of density flux
	586	!-- in the bottom and top surface layer
	587	DO j = nys, nyn
[2232]	588	DO k = nzb+1, nzt
[2126]	589	tend(k,j,i) = tend(k,j,i) + &
	590	kh(k,j,i) * g / prho(k,j,i) * &
	591	( prho(k+1,j,i) - prho(k-1,j,i) ) * &
[2232]	592	dd2zu(k) * &
	593	MERGE( 1.0_wp, 0.0_wp, &
	594	BTEST( wall_flags_0(k,j,i), 0 ) &
	595	)
[940]	596	ENDDO
[97]	597	ENDDO
	598
[940]	599	ELSE
[97]	600
[940]	601	DO j = nys, nyn
[2232]	602	DO k = nzb+1, nzt
	603	!
	604	!-- Flag 9 is used to mask top fluxes, flag 30 to mask
	605	!-- surface fluxes
	606	tend(k,j,i) = tend(k,j,i) - &
	607	kh(k,j,i) * g / pt(k,j,i) * &
	608	( pt(k+1,j,i) - pt(k-1,j,i) ) * &
	609	dd2zu(k) * &
	610	MERGE( 1.0_wp, 0.0_wp, &
	611	BTEST( wall_flags_0(k,j,i), 30 ) &
	612	) * &
	613	MERGE( 1.0_wp, 0.0_wp, &
	614	BTEST( wall_flags_0(k,j,i), 9 ) &
	615	)
[940]	616	ENDDO
	617
	618	IF ( use_surface_fluxes ) THEN
[2232]	619	!
	620	!-- Default surfaces, up- and downwrd-facing
	621	DO l = 0, 1
	622	surf_s = surf_def_h(l)%start_index(j,i)
	623	surf_e = surf_def_h(l)%end_index(j,i)
	624	DO m = surf_s, surf_e
	625	k = surf_def_h(l)%k(m)
	626	tend(k,j,i) = tend(k,j,i) + g / pt_reference &
[2329]	627	* drho_air_zw(k-1) &
[2232]	628	* surf_def_h(l)%shf(m)
	629	ENDDO
	630	ENDDO
	631	!
	632	!-- Natural surfaces
	633	surf_s = surf_lsm_h%start_index(j,i)
	634	surf_e = surf_lsm_h%end_index(j,i)
	635	DO m = surf_s, surf_e
	636	k = surf_lsm_h%k(m)
	637	tend(k,j,i) = tend(k,j,i) + g / pt_reference &
[2329]	638	* drho_air_zw(k-1) &
[2232]	639	* surf_lsm_h%shf(m)
	640	ENDDO
	641	!
	642	!-- Urban surfaces
	643	surf_s = surf_usm_h%start_index(j,i)
	644	surf_e = surf_usm_h%end_index(j,i)
	645	DO m = surf_s, surf_e
	646	k = surf_usm_h%k(m)
	647	tend(k,j,i) = tend(k,j,i) + g / pt_reference &
[2329]	648	* drho_air_zw(k-1) &
[2232]	649	* surf_usm_h%shf(m)
	650	ENDDO
[940]	651	ENDIF
	652
	653	IF ( use_top_fluxes ) THEN
[2232]	654	surf_s = surf_def_h(2)%start_index(j,i)
	655	surf_e = surf_def_h(2)%end_index(j,i)
	656	DO m = surf_s, surf_e
	657	k = surf_def_h(2)%k(m)
[2329]	658	tend(k,j,i) = tend(k,j,i) + g / pt_reference &
	659	* drho_air_zw(k-1) &
	660	* surf_def_h(2)%shf(m)
[2232]	661	ENDDO
[940]	662	ENDIF
[97]	663	ENDDO
	664
[940]	665	ENDIF
[97]	666
	667	ENDIF
[1]	668
[940]	669	ELSE
[57]	670
[940]	671	DO j = nys, nyn
[1]	672
[2232]	673	DO k = nzb+1, nzt
	674	!
	675	!-- Flag 9 is used to mask top fluxes, flag 30 to mask
	676	!-- surface fluxes
[1007]	677	IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
[1342]	678	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
	679	k2 = 0.61_wp * pt(k,j,i)
[1007]	680	tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * &
	681	g / vpt(k,j,i) * &
	682	( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + &
	683	k2 * ( q(k+1,j,i) - q(k-1,j,i) ) &
[2232]	684	) * dd2zu(k) * &
	685	MERGE( 1.0_wp, 0.0_wp, &
	686	BTEST( wall_flags_0(k,j,i), 30 ) &
	687	) * &
	688	MERGE( 1.0_wp, 0.0_wp, &
	689	BTEST( wall_flags_0(k,j,i), 9 ) &
	690	)
[1007]	691	ELSE IF ( cloud_physics ) THEN
[1342]	692	IF ( ql(k,j,i) == 0.0_wp ) THEN
	693	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
	694	k2 = 0.61_wp * pt(k,j,i)
[940]	695	ELSE
	696	theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
	697	temp = theta * t_d_pt(k)
[1342]	698	k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * &
[2232]	699	( q(k,j,i) - ql(k,j,i) ) * &
	700	( 1.0_wp + 0.622_wp * l_d_r / temp ) ) / &
	701	( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * &
[940]	702	( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
[1342]	703	k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
[940]	704	ENDIF
[1007]	705	tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * &
	706	g / vpt(k,j,i) * &
[940]	707	( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + &
	708	k2 * ( q(k+1,j,i) - q(k-1,j,i) ) &
[2232]	709	) * dd2zu(k) * &
	710	MERGE( 1.0_wp, 0.0_wp, &
	711	BTEST( wall_flags_0(k,j,i), 30 ) &
	712	) * &
	713	MERGE( 1.0_wp, 0.0_wp, &
	714	BTEST( wall_flags_0(k,j,i), 9 ) &
	715	)
[1007]	716	ELSE IF ( cloud_droplets ) THEN
[1342]	717	k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
	718	k2 = 0.61_wp * pt(k,j,i)
[2232]	719	tend(k,j,i) = tend(k,j,i) - &
	720	kh(k,j,i) * g / vpt(k,j,i) * &
	721	( k1 * ( pt(k+1,j,i)- pt(k-1,j,i) ) + &
	722	k2 * ( q(k+1,j,i) - q(k-1,j,i) ) - &
	723	pt(k,j,i) * ( ql(k+1,j,i) - &
	724	ql(k-1,j,i) ) ) * dd2zu(k) * &
	725	MERGE( 1.0_wp, 0.0_wp, &
	726	BTEST( wall_flags_0(k,j,i), 30 ) &
	727	) * &
	728	MERGE( 1.0_wp, 0.0_wp, &
	729	BTEST( wall_flags_0(k,j,i), 9 ) &
	730	)
[1007]	731	ENDIF
	732
[940]	733	ENDDO
	734
[1]	735	ENDDO
	736
[940]	737	IF ( use_surface_fluxes ) THEN
[1]	738
[940]	739	DO j = nys, nyn
[2232]	740	!
	741	!-- Treat horizontal default surfaces, up- and downward-facing
	742	DO l = 0, 1
	743	surf_s = surf_def_h(l)%start_index(j,i)
	744	surf_e = surf_def_h(l)%end_index(j,i)
	745	DO m = surf_s, surf_e
	746	k = surf_def_h(l)%k(m)
[1]	747
[2232]	748	IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
	749	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
	750	k2 = 0.61_wp * pt(k,j,i)
	751	ELSE IF ( cloud_physics ) THEN
	752	IF ( ql(k,j,i) == 0.0_wp ) THEN
	753	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
	754	k2 = 0.61_wp * pt(k,j,i)
	755	ELSE
	756	theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
	757	temp = theta * t_d_pt(k)
	758	k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * &
	759	( q(k,j,i) - ql(k,j,i) ) * &
	760	( 1.0_wp + 0.622_wp * l_d_r / temp ) ) / &
	761	( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * &
	762	( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
	763	k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
	764	ENDIF
	765	ELSE IF ( cloud_droplets ) THEN
	766	k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
	767	k2 = 0.61_wp * pt(k,j,i)
	768	ENDIF
[1]	769
[2232]	770	tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * &
	771	( k1 * surf_def_h(l)%shf(m) + &
[2329]	772	k2 * surf_def_h(l)%qsws(m) &
	773	) * drho_air_zw(k-1)
[2232]	774	ENDDO
	775	ENDDO
	776	!
	777	!-- Treat horizontal natural surfaces
	778	surf_s = surf_lsm_h%start_index(j,i)
	779	surf_e = surf_lsm_h%end_index(j,i)
	780	DO m = surf_s, surf_e
	781	k = surf_lsm_h%k(m)
	782
	783	IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
[1342]	784	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
	785	k2 = 0.61_wp * pt(k,j,i)
[2232]	786	ELSE IF ( cloud_physics ) THEN
	787	IF ( ql(k,j,i) == 0.0_wp ) THEN
	788	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
	789	k2 = 0.61_wp * pt(k,j,i)
	790	ELSE
	791	theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
	792	temp = theta * t_d_pt(k)
	793	k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * &
	794	( q(k,j,i) - ql(k,j,i) ) * &
	795	( 1.0_wp + 0.622_wp * l_d_r / temp ) ) / &
	796	( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * &
	797	( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
	798	k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
	799	ENDIF
	800	ELSE IF ( cloud_droplets ) THEN
	801	k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
	802	k2 = 0.61_wp * pt(k,j,i)
[940]	803	ENDIF
[1]	804
[2232]	805	tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * &
	806	( k1 * surf_lsm_h%shf(m) + &
[2329]	807	k2 * surf_lsm_h%qsws(m) &
	808	) * drho_air_zw(k-1)
[2232]	809	ENDDO
	810	!
	811	!-- Treat horizontal urban surfaces
	812	surf_s = surf_usm_h%start_index(j,i)
	813	surf_e = surf_usm_h%end_index(j,i)
	814	DO m = surf_s, surf_e
	815	k = surf_lsm_h%k(m)
	816
	817	IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
	818	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
	819	k2 = 0.61_wp * pt(k,j,i)
	820	ELSE IF ( cloud_physics ) THEN
	821	IF ( ql(k,j,i) == 0.0_wp ) THEN
	822	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
	823	k2 = 0.61_wp * pt(k,j,i)
	824	ELSE
	825	theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
	826	temp = theta * t_d_pt(k)
	827	k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * &
	828	( q(k,j,i) - ql(k,j,i) ) * &
	829	( 1.0_wp + 0.622_wp * l_d_r / temp ) ) / &
	830	( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * &
	831	( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
	832	k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
	833	ENDIF
	834	ELSE IF ( cloud_droplets ) THEN
	835	k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
	836	k2 = 0.61_wp * pt(k,j,i)
	837	ENDIF
	838
	839	tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * &
	840	( k1 * surf_usm_h%shf(m) + &
[2329]	841	k2 * surf_usm_h%qsws(m) &
	842	) * drho_air_zw(k-1)
[2232]	843	ENDDO
	844
[940]	845	ENDDO
[1]	846
[940]	847	ENDIF
[1]	848
[940]	849	IF ( use_top_fluxes ) THEN
[19]	850
[940]	851	DO j = nys, nyn
[19]	852
[2232]	853	surf_s = surf_def_h(2)%start_index(j,i)
	854	surf_e = surf_def_h(2)%end_index(j,i)
	855	DO m = surf_s, surf_e
	856	k = surf_def_h(2)%k(m)
[19]	857
[2232]	858	IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
[1342]	859	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
	860	k2 = 0.61_wp * pt(k,j,i)
[2232]	861	ELSE IF ( cloud_physics ) THEN
	862	IF ( ql(k,j,i) == 0.0_wp ) THEN
	863	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
	864	k2 = 0.61_wp * pt(k,j,i)
	865	ELSE
	866	theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
	867	temp = theta * t_d_pt(k)
	868	k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * &
[1353]	869	( q(k,j,i) - ql(k,j,i) ) * &
	870	( 1.0_wp + 0.622_wp * l_d_r / temp ) ) / &
	871	( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * &
[940]	872	( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
[2232]	873	k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
	874	ENDIF
	875	ELSE IF ( cloud_droplets ) THEN
	876	k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
	877	k2 = 0.61_wp * pt(k,j,i)
[940]	878	ENDIF
[19]	879
[2232]	880	tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * &
	881	( k1 * surf_def_h(2)%shf(m) + &
[2329]	882	k2 * surf_def_h(2)%qsws(m) &
	883	) * drho_air_zw(k-1)
[2232]	884
	885	ENDDO
	886
[940]	887	ENDDO
[19]	888
[940]	889	ENDIF
	890
[19]	891	ENDIF
	892
[1]	893	ENDIF
	894
	895	ENDDO
	896
	897	END SUBROUTINE production_e
	898
	899
	900	!------------------------------------------------------------------------------!
[1682]	901	! Description:
	902	! ------------
	903	!> Call for grid point i,j
[1]	904	!------------------------------------------------------------------------------!
	905	SUBROUTINE production_e_ij( i, j )
	906
[1320]	907	USE arrays_3d, &
[2329]	908	ONLY: ddzw, dd2zu, drho_air_zw, kh, km, prho, pt, q, ql, tend, u, &
	909	v, vpt, w
[449]	910
[1320]	911	USE cloud_parameters, &
	912	ONLY: l_d_cp, l_d_r, pt_d_t, t_d_pt
	913
	914	USE control_parameters, &
[1691]	915	ONLY: cloud_droplets, cloud_physics, constant_flux_layer, g, &
	916	humidity, kappa, neutral, ocean, pt_reference, &
	917	rho_reference, use_single_reference_value, &
	918	use_surface_fluxes, use_top_fluxes
[1320]	919
	920	USE grid_variables, &
[2232]	921	ONLY: ddx, dx, ddy, dy
[1320]	922
	923	USE indices, &
[2232]	924	ONLY: nxl, nxr, nys, nyn, nzb, nzb, nzt, wall_flags_0
[1320]	925
[2232]	926	USE surface_mod, &
	927	ONLY : surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, &
	928	surf_usm_v
	929
[1]	930	IMPLICIT NONE
	931
[2232]	932	INTEGER(iwp) :: i !< running index x-direction
	933	INTEGER(iwp) :: j !< running index y-direction
	934	INTEGER(iwp) :: k !< running index z-direction
	935	INTEGER(iwp) :: l !< running index for different surface type orientation
	936	INTEGER(iwp) :: m !< running index surface elements
	937	INTEGER(iwp) :: surf_e !< end index of surface elements at given i-j position
	938	INTEGER(iwp) :: surf_s !< start index of surface elements at given i-j position
[1]	939
[1682]	940	REAL(wp) :: def !<
[2232]	941	REAL(wp) :: flag !< flag to mask topography
[1682]	942	REAL(wp) :: k1 !<
	943	REAL(wp) :: k2 !<
[2232]	944	REAL(wp) :: km_neutral !< diffusion coefficient assuming neutral conditions - used to compute shear production at surfaces
[1682]	945	REAL(wp) :: theta !<
	946	REAL(wp) :: temp !<
[2232]	947	REAL(wp) :: sign_dir !< sign of wall-tke flux, depending on wall orientation
	948	REAL(wp) :: usvs !< momentum flux u"v"
	949	REAL(wp) :: vsus !< momentum flux v"u"
	950	REAL(wp) :: wsus !< momentum flux w"u"
	951	REAL(wp) :: wsvs !< momentum flux w"v"
[1]	952
[53]	953
[2232]	954	REAL(wp), DIMENSION(nzb+1:nzt) :: dudx !< Gradient of u-component in x-direction
	955	REAL(wp), DIMENSION(nzb+1:nzt) :: dudy !< Gradient of u-component in y-direction
	956	REAL(wp), DIMENSION(nzb+1:nzt) :: dudz !< Gradient of u-component in z-direction
	957	REAL(wp), DIMENSION(nzb+1:nzt) :: dvdx !< Gradient of v-component in x-direction
	958	REAL(wp), DIMENSION(nzb+1:nzt) :: dvdy !< Gradient of v-component in y-direction
	959	REAL(wp), DIMENSION(nzb+1:nzt) :: dvdz !< Gradient of v-component in z-direction
	960	REAL(wp), DIMENSION(nzb+1:nzt) :: dwdx !< Gradient of w-component in x-direction
	961	REAL(wp), DIMENSION(nzb+1:nzt) :: dwdy !< Gradient of w-component in y-direction
	962	REAL(wp), DIMENSION(nzb+1:nzt) :: dwdz !< Gradient of w-component in z-direction
[1]	963
	964
[2232]	965	IF ( constant_flux_layer ) THEN
	966	!
	967	!-- Calculate TKE production by shear. Calculate gradients at all grid
	968	!-- points first, gradients at surface-bounded grid points will be
	969	!-- overwritten further below.
	970	DO k = nzb+1, nzt
[1]	971
[2232]	972	dudx(k) = ( u(k,j,i+1) - u(k,j,i) ) * ddx
	973	dudy(k) = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
	974	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
	975	dudz(k) = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - &
	976	u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
[1]	977
[2232]	978	dvdx(k) = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
	979	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
	980	dvdy(k) = ( v(k,j+1,i) - v(k,j,i) ) * ddy
	981	dvdz(k) = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - &
	982	v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
[1]	983
[2232]	984	dwdx(k) = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
	985	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
	986	dwdy(k) = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
	987	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
	988	dwdz(k) = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
[1]	989
[2232]	990	ENDDO
	991	!
	992	!-- Compute gradients at north- and south-facing surfaces.
	993	!-- Note, no vertical natural surfaces so far.
	994	DO l = 0, 1
	995	!
	996	!-- Default surfaces
	997	surf_s = surf_def_v(l)%start_index(j,i)
	998	surf_e = surf_def_v(l)%end_index(j,i)
	999	DO m = surf_s, surf_e
	1000	k = surf_def_v(l)%k(m)
	1001	usvs = surf_def_v(l)%mom_flux_tke(0,m)
	1002	wsvs = surf_def_v(l)%mom_flux_tke(1,m)
[1007]	1003
[2232]	1004	km_neutral = kappa * ( usvs2 + wsvs2 )**0.25_wp &
	1005	* 0.5_wp * dy
[1]	1006	!
[2232]	1007	!-- -1.0 for right-facing wall, 1.0 for left-facing wall
	1008	sign_dir = MERGE( 1.0_wp, -1.0_wp, &
	1009	BTEST( wall_flags_0(k,j-1,i), 0 ) )
	1010	dudy(k) = sign_dir * usvs / ( km_neutral + 1E-10_wp )
	1011	dwdy(k) = sign_dir * wsvs / ( km_neutral + 1E-10_wp )
	1012	ENDDO
	1013	!
	1014	!-- Natural surfaces
	1015	surf_s = surf_lsm_v(l)%start_index(j,i)
	1016	surf_e = surf_lsm_v(l)%end_index(j,i)
	1017	DO m = surf_s, surf_e
	1018	k = surf_lsm_v(l)%k(m)
	1019	usvs = surf_lsm_v(l)%mom_flux_tke(0,m)
	1020	wsvs = surf_lsm_v(l)%mom_flux_tke(1,m)
[1]	1021
[2232]	1022	km_neutral = kappa * ( usvs2 + wsvs2 )**0.25_wp &
	1023	* 0.5_wp * dy
[1007]	1024	!
[2232]	1025	!-- -1.0 for right-facing wall, 1.0 for left-facing wall
	1026	sign_dir = MERGE( 1.0_wp, -1.0_wp, &
	1027	BTEST( wall_flags_0(k,j-1,i), 0 ) )
	1028	dudy(k) = sign_dir * usvs / ( km_neutral + 1E-10_wp )
	1029	dwdy(k) = sign_dir * wsvs / ( km_neutral + 1E-10_wp )
	1030	ENDDO
	1031	!
	1032	!-- Urban surfaces
	1033	surf_s = surf_usm_v(l)%start_index(j,i)
	1034	surf_e = surf_usm_v(l)%end_index(j,i)
	1035	DO m = surf_s, surf_e
	1036	k = surf_usm_v(l)%k(m)
	1037	usvs = surf_usm_v(l)%mom_flux_tke(0,m)
	1038	wsvs = surf_usm_v(l)%mom_flux_tke(1,m)
[1]	1039
[2232]	1040	km_neutral = kappa * ( usvs2 + wsvs2 )**0.25_wp &
	1041	* 0.5_wp * dy
[1007]	1042	!
[2232]	1043	!-- -1.0 for right-facing wall, 1.0 for left-facing wall
	1044	sign_dir = MERGE( 1.0_wp, -1.0_wp, &
	1045	BTEST( wall_flags_0(k,j-1,i), 0 ) )
	1046	dudy(k) = sign_dir * usvs / ( km_neutral + 1E-10_wp )
	1047	dwdy(k) = sign_dir * wsvs / ( km_neutral + 1E-10_wp )
	1048	ENDDO
	1049	ENDDO
	1050	!
	1051	!-- Compute gradients at east- and west-facing walls
	1052	DO l = 2, 3
	1053	!
	1054	!-- Default surfaces
	1055	surf_s = surf_def_v(l)%start_index(j,i)
	1056	surf_e = surf_def_v(l)%end_index(j,i)
	1057	DO m = surf_s, surf_e
	1058	k = surf_def_v(l)%k(m)
	1059	vsus = surf_def_v(l)%mom_flux_tke(0,m)
	1060	wsus = surf_def_v(l)%mom_flux_tke(1,m)
[1]	1061
[2232]	1062	km_neutral = kappa * ( vsus2 + wsus2 )**0.25_wp &
	1063	* 0.5_wp * dx
[1]	1064	!
[2232]	1065	!-- -1.0 for right-facing wall, 1.0 for left-facing wall
	1066	sign_dir = MERGE( 1.0_wp, -1.0_wp, &
	1067	BTEST( wall_flags_0(k,j,i-1), 0 ) )
	1068	dvdx(k) = sign_dir * vsus / ( km_neutral + 1E-10_wp )
	1069	dwdx(k) = sign_dir * wsus / ( km_neutral + 1E-10_wp )
	1070	ENDDO
	1071	!
	1072	!-- Natural surfaces
	1073	surf_s = surf_lsm_v(l)%start_index(j,i)
	1074	surf_e = surf_lsm_v(l)%end_index(j,i)
	1075	DO m = surf_s, surf_e
	1076	k = surf_lsm_v(l)%k(m)
	1077	vsus = surf_lsm_v(l)%mom_flux_tke(0,m)
	1078	wsus = surf_lsm_v(l)%mom_flux_tke(1,m)
[1]	1079
[2232]	1080	km_neutral = kappa * ( vsus2 + wsus2 )**0.25_wp &
	1081	* 0.5_wp * dx
[1007]	1082	!
[2232]	1083	!-- -1.0 for right-facing wall, 1.0 for left-facing wall
	1084	sign_dir = MERGE( 1.0_wp, -1.0_wp, &
	1085	BTEST( wall_flags_0(k,j,i-1), 0 ) )
	1086	dvdx(k) = sign_dir * vsus / ( km_neutral + 1E-10_wp )
	1087	dwdx(k) = sign_dir * wsus / ( km_neutral + 1E-10_wp )
	1088	ENDDO
	1089	!
	1090	!-- Urban surfaces
	1091	surf_s = surf_usm_v(l)%start_index(j,i)
	1092	surf_e = surf_usm_v(l)%end_index(j,i)
	1093	DO m = surf_s, surf_e
	1094	k = surf_usm_v(l)%k(m)
	1095	vsus = surf_usm_v(l)%mom_flux_tke(0,m)
	1096	wsus = surf_usm_v(l)%mom_flux_tke(1,m)
[1]	1097
[2232]	1098	km_neutral = kappa * ( vsus2 + wsus2 )**0.25_wp &
	1099	* 0.5_wp * dx
[1007]	1100	!
[2232]	1101	!-- -1.0 for right-facing wall, 1.0 for left-facing wall
	1102	sign_dir = MERGE( 1.0_wp, -1.0_wp, &
	1103	BTEST( wall_flags_0(k,j,i-1), 0 ) )
	1104	dvdx(k) = sign_dir * vsus / ( km_neutral + 1E-10_wp )
	1105	dwdx(k) = sign_dir * wsus / ( km_neutral + 1E-10_wp )
	1106	ENDDO
	1107	ENDDO
	1108	!
	1109	!-- Compute gradients at upward-facing walls, first for
	1110	!-- non-natural default surfaces
	1111	surf_s = surf_def_h(0)%start_index(j,i)
	1112	surf_e = surf_def_h(0)%end_index(j,i)
	1113	DO m = surf_s, surf_e
	1114	k = surf_def_h(0)%k(m)
	1115	!
	1116	!-- Please note, actually, an interpolation of u_0 and v_0
	1117	!-- onto the grid center would be required. However, this
	1118	!-- would require several data transfers between 2D-grid and
	1119	!-- wall type. The effect of this missing interpolation is
	1120	!-- negligible. (See also production_e_init).
	1121	dudz(k) = ( u(k+1,j,i) - surf_def_h(0)%u_0(m) ) * dd2zu(k)
	1122	dvdz(k) = ( v(k+1,j,i) - surf_def_h(0)%v_0(m) ) * dd2zu(k)
[1]	1123
[2232]	1124	ENDDO
[1]	1125	!
[2232]	1126	!-- Natural surfaces
	1127	surf_s = surf_lsm_h%start_index(j,i)
	1128	surf_e = surf_lsm_h%end_index(j,i)
	1129	DO m = surf_s, surf_e
	1130	k = surf_lsm_h%k(m)
	1131	!
	1132	!-- Please note, actually, an interpolation of u_0 and v_0
	1133	!-- onto the grid center would be required. However, this
	1134	!-- would require several data transfers between 2D-grid and
	1135	!-- wall type. The effect of this missing interpolation is
	1136	!-- negligible. (See also production_e_init).
	1137	dudz(k) = ( u(k+1,j,i) - surf_lsm_h%u_0(m) ) * dd2zu(k)
	1138	dvdz(k) = ( v(k+1,j,i) - surf_lsm_h%v_0(m) ) * dd2zu(k)
	1139	ENDDO
	1140	!
	1141	!-- Urban surfaces
	1142	surf_s = surf_usm_h%start_index(j,i)
	1143	surf_e = surf_usm_h%end_index(j,i)
	1144	DO m = surf_s, surf_e
	1145	k = surf_usm_h%k(m)
	1146	!
	1147	!-- Please note, actually, an interpolation of u_0 and v_0
	1148	!-- onto the grid center would be required. However, this
	1149	!-- would require several data transfers between 2D-grid and
	1150	!-- wall type. The effect of this missing interpolation is
	1151	!-- negligible. (See also production_e_init).
	1152	dudz(k) = ( u(k+1,j,i) - surf_usm_h%u_0(m) ) * dd2zu(k)
	1153	dvdz(k) = ( v(k+1,j,i) - surf_usm_h%v_0(m) ) * dd2zu(k)
	1154	ENDDO
	1155	!
	1156	!-- Compute gradients at downward-facing walls, only for
	1157	!-- non-natural default surfaces
	1158	surf_s = surf_def_h(1)%start_index(j,i)
	1159	surf_e = surf_def_h(1)%end_index(j,i)
	1160	DO m = surf_s, surf_e
	1161	k = surf_def_h(1)%k(m)
	1162	!
	1163	!-- Please note, actually, an interpolation of u_0 and v_0
	1164	!-- onto the grid center would be required. However, this
	1165	!-- would require several data transfers between 2D-grid and
	1166	!-- wall type. The effect of this missing interpolation is
	1167	!-- negligible. (See also production_e_init).
	1168	dudz(k) = ( surf_def_h(1)%u_0(m) - u(k-1,j,i) ) * dd2zu(k)
	1169	dvdz(k) = ( surf_def_h(1)%v_0(m) - v(k-1,j,i) ) * dd2zu(k)
[1]	1170
[2232]	1171	ENDDO
[1]	1172
[2232]	1173	DO k = nzb+1, nzt
[1]	1174
[2232]	1175	def = 2.0_wp * ( dudx(k)2 + dvdy(k)2 + dwdz(k)**2 ) + &
	1176	dudy(k)2 + dvdx(k)2 + dwdx(k)**2 + &
	1177	dwdy(k)2 + dudz(k)2 + dvdz(k)**2 + &
	1178	2.0_wp * ( dvdx(k)dudy(k) + dwdx(k)dudz(k) + dwdy(k)*dvdz(k) )
[1]	1179
[1342]	1180	IF ( def < 0.0_wp ) def = 0.0_wp
[1]	1181
[2232]	1182	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
[1]	1183
[2232]	1184	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def * flag
[1]	1185
[2232]	1186	ENDDO
	1187
[2478]	1188	ELSE
[1]	1189	!
[2232]	1190	!-- Calculate TKE production by shear. Here, no additional
	1191	!-- wall-bounded code is considered.
	1192	!-- Why?
	1193	DO k = nzb+1, nzt
[1]	1194
[2232]	1195	dudx(k) = ( u(k,j,i+1) - u(k,j,i) ) * ddx
	1196	dudy(k) = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
	1197	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
	1198	dudz(k) = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - &
	1199	u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
[1]	1200
[2232]	1201	dvdx(k) = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
	1202	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
	1203	dvdy(k) = ( v(k,j+1,i) - v(k,j,i) ) * ddy
	1204	dvdz(k) = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - &
	1205	v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
[1]	1206
[2232]	1207	dwdx(k) = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
	1208	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
	1209	dwdy(k) = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
	1210	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
	1211	dwdz(k) = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
[1]	1212
[2232]	1213	def = 2.0_wp * ( dudx(k)2 + dvdy(k)2 + dwdz(k)**2 ) + &
	1214	dudy(k)2 + dvdx(k)2 + dwdx(k)**2 + &
	1215	dwdy(k)2 + dudz(k)2 + dvdz(k)**2 + &
	1216	2.0_wp * ( dvdx(k)dudy(k) + dwdx(k)dudz(k) + dwdy(k)*dvdz(k) )
[1]	1217
[1342]	1218	IF ( def < 0.0_wp ) def = 0.0_wp
[1]	1219
[2232]	1220	flag = MERGE( 1.0_wp, 0.0_wp, &
	1221	BTEST( wall_flags_0(k,j,i), 29 ) )
	1222	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def * flag
[1]	1223
[2232]	1224	ENDDO
[1]	1225
	1226	ENDIF
	1227
	1228	!
[940]	1229	!-- If required, calculate TKE production by buoyancy
	1230	IF ( .NOT. neutral ) THEN
[1]	1231
[940]	1232	IF ( .NOT. humidity ) THEN
[19]	1233
[1179]	1234	IF ( use_single_reference_value ) THEN
[940]	1235
	1236	IF ( ocean ) THEN
[97]	1237	!
[940]	1238	!-- So far in the ocean no special treatment of density flux in
	1239	!-- the bottom and top surface layer
[2232]	1240	DO k = nzb+1, nzt
[2478]	1241
[2232]	1242	tend(k,j,i) = tend(k,j,i) + &
	1243	kh(k,j,i) * g / rho_reference * &
	1244	( prho(k+1,j,i) - prho(k-1,j,i) ) * &
	1245	dd2zu(k) * &
	1246	MERGE( 1.0_wp, 0.0_wp, &
	1247	BTEST( wall_flags_0(k,j,i), 0 ) &
[2478]	1248	)
[940]	1249	ENDDO
[97]	1250
[940]	1251	ELSE
[97]	1252
[2232]	1253	DO k = nzb+1, nzt
	1254	!
	1255	!-- Flag 9 is used to mask top fluxes, flag 30 to mask
	1256	!-- surface fluxes
	1257	tend(k,j,i) = tend(k,j,i) - &
	1258	kh(k,j,i) * g / pt_reference * &
	1259	( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k) * &
	1260	MERGE( 1.0_wp, 0.0_wp, &
	1261	BTEST( wall_flags_0(k,j,i), 30 ) &
	1262	) * &
	1263	MERGE( 1.0_wp, 0.0_wp, &
	1264	BTEST( wall_flags_0(k,j,i), 9 ) &
	1265	)
	1266
[940]	1267	ENDDO
[1]	1268
[940]	1269	IF ( use_surface_fluxes ) THEN
[2232]	1270	!
	1271	!-- Default surfaces, up- and downward-facing
	1272	DO l = 0, 1
	1273	surf_s = surf_def_h(l)%start_index(j,i)
	1274	surf_e = surf_def_h(l)%end_index(j,i)
	1275	DO m = surf_s, surf_e
	1276	k = surf_def_h(l)%k(m)
	1277	tend(k,j,i) = tend(k,j,i) + g / pt_reference * &
[2329]	1278	drho_air_zw(k-1) * &
	1279	surf_def_h(l)%shf(m)
[2232]	1280	ENDDO
	1281	ENDDO
	1282	!
	1283	!-- Natural surfaces
	1284	surf_s = surf_lsm_h%start_index(j,i)
	1285	surf_e = surf_lsm_h%end_index(j,i)
	1286	DO m = surf_s, surf_e
	1287	k = surf_lsm_h%k(m)
	1288	tend(k,j,i) = tend(k,j,i) + g / pt_reference * &
[2329]	1289	drho_air_zw(k-1) * &
[2232]	1290	surf_lsm_h%shf(m)
	1291	ENDDO
	1292	!
	1293	!-- Urban surfaces
	1294	surf_s = surf_usm_h%start_index(j,i)
	1295	surf_e = surf_usm_h%end_index(j,i)
	1296	DO m = surf_s, surf_e
	1297	k = surf_usm_h%k(m)
	1298	tend(k,j,i) = tend(k,j,i) + g / pt_reference * &
[2329]	1299	drho_air_zw(k-1) * &
[2232]	1300	surf_usm_h%shf(m)
	1301	ENDDO
[940]	1302	ENDIF
[19]	1303
[940]	1304	IF ( use_top_fluxes ) THEN
[2232]	1305	surf_s = surf_def_h(2)%start_index(j,i)
	1306	surf_e = surf_def_h(2)%end_index(j,i)
	1307	DO m = surf_s, surf_e
	1308	k = surf_def_h(2)%k(m)
	1309	tend(k,j,i) = tend(k,j,i) + g / pt_reference * &
[2329]	1310	drho_air_zw(k-1) * &
[2232]	1311	surf_def_h(2)%shf(m)
	1312	ENDDO
[940]	1313	ENDIF
	1314
[97]	1315	ENDIF
	1316
[940]	1317	ELSE
[57]	1318
[940]	1319	IF ( ocean ) THEN
[97]	1320	!
[940]	1321	!-- So far in the ocean no special treatment of density flux in
	1322	!-- the bottom and top surface layer
[2232]	1323	DO k = nzb+1, nzt
[2126]	1324	tend(k,j,i) = tend(k,j,i) + &
[2232]	1325	kh(k,j,i) * g / prho(k,j,i) * &
	1326	( prho(k+1,j,i) - prho(k-1,j,i) ) * &
	1327	dd2zu(k) * &
	1328	MERGE( 1.0_wp, 0.0_wp, &
	1329	BTEST( wall_flags_0(k,j,i), 0 ) &
[2478]	1330	)
[940]	1331	ENDDO
[97]	1332
[940]	1333	ELSE
[97]	1334
[2232]	1335	DO k = nzb+1, nzt
	1336	!
	1337	!-- Flag 9 is used to mask top fluxes, flag 30 to mask
	1338	!-- surface fluxes
	1339	tend(k,j,i) = tend(k,j,i) - &
	1340	kh(k,j,i) * g / pt(k,j,i) * &
	1341	( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k) * &
	1342	MERGE( 1.0_wp, 0.0_wp, &
	1343	BTEST( wall_flags_0(k,j,i), 30 ) &
	1344	) * &
	1345	MERGE( 1.0_wp, 0.0_wp, &
	1346	BTEST( wall_flags_0(k,j,i), 9 ) &
	1347	)
[940]	1348	ENDDO
[57]	1349
[940]	1350	IF ( use_surface_fluxes ) THEN
[2232]	1351	!
	1352	!-- Default surfaces, up- and downward-facing
	1353	DO l = 0, 1
	1354	surf_s = surf_def_h(l)%start_index(j,i)
	1355	surf_e = surf_def_h(l)%end_index(j,i)
	1356	DO m = surf_s, surf_e
	1357	k = surf_def_h(l)%k(m)
	1358	tend(k,j,i) = tend(k,j,i) + g / pt_reference &
[2329]	1359	* drho_air_zw(k-1) &
	1360	* surf_def_h(l)%shf(m)
[2232]	1361	ENDDO
	1362	ENDDO
	1363	!
	1364	!-- Natural surfaces
	1365	surf_s = surf_lsm_h%start_index(j,i)
	1366	surf_e = surf_lsm_h%end_index(j,i)
	1367	DO m = surf_s, surf_e
	1368	k = surf_lsm_h%k(m)
	1369	tend(k,j,i) = tend(k,j,i) + g / pt_reference &
[2329]	1370	* drho_air_zw(k-1) &
	1371	* surf_lsm_h%shf(m)
[2232]	1372	ENDDO
	1373	!
	1374	!-- Urban surfaces
	1375	surf_s = surf_usm_h%start_index(j,i)
	1376	surf_e = surf_usm_h%end_index(j,i)
	1377	DO m = surf_s, surf_e
	1378	k = surf_usm_h%k(m)
	1379	tend(k,j,i) = tend(k,j,i) + g / pt_reference &
[2329]	1380	* drho_air_zw(k-1) &
	1381	* surf_usm_h%shf(m)
[2232]	1382	ENDDO
[940]	1383	ENDIF
[57]	1384
[940]	1385	IF ( use_top_fluxes ) THEN
[2232]	1386	surf_s = surf_def_h(2)%start_index(j,i)
	1387	surf_e = surf_def_h(2)%end_index(j,i)
	1388	DO m = surf_s, surf_e
	1389	k = surf_def_h(2)%k(m)
	1390	tend(k,j,i) = tend(k,j,i) + g / pt_reference * &
[2329]	1391	drho_air_zw(k-1) * &
	1392	surf_def_h(2)%shf(m)
[2232]	1393	ENDDO
[940]	1394	ENDIF
	1395
[97]	1396	ENDIF
	1397
[57]	1398	ENDIF
	1399
[940]	1400	ELSE
[57]	1401
[2232]	1402	DO k = nzb+1, nzt
	1403	!
	1404	!-- Flag 9 is used to mask top fluxes, flag 30 to mask surface fluxes
[1007]	1405	IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
[1342]	1406	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
	1407	k2 = 0.61_wp * pt(k,j,i)
[2232]	1408	tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) * &
	1409	( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + &
	1410	k2 * ( q(k+1,j,i) - q(k-1,j,i) ) &
	1411	) * dd2zu(k) * &
	1412	MERGE( 1.0_wp, 0.0_wp, &
	1413	BTEST( wall_flags_0(k,j,i), 30 ) &
	1414	) * &
	1415	MERGE( 1.0_wp, 0.0_wp, &
	1416	BTEST( wall_flags_0(k,j,i), 9 ) &
	1417	)
	1418
[1007]	1419	ELSE IF ( cloud_physics ) THEN
[1342]	1420	IF ( ql(k,j,i) == 0.0_wp ) THEN
	1421	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
	1422	k2 = 0.61_wp * pt(k,j,i)
[940]	1423	ELSE
	1424	theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
	1425	temp = theta * t_d_pt(k)
[2232]	1426	k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * &
	1427	( q(k,j,i) - ql(k,j,i) ) * &
	1428	( 1.0_wp + 0.622_wp * l_d_r / temp ) ) / &
	1429	( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * &
[940]	1430	( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
[1342]	1431	k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
[940]	1432	ENDIF
[2232]	1433	tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) * &
	1434	( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + &
	1435	k2 * ( q(k+1,j,i) - q(k-1,j,i) ) &
	1436	) * dd2zu(k) * &
	1437	MERGE( 1.0_wp, 0.0_wp, &
	1438	BTEST( wall_flags_0(k,j,i), 30 ) &
	1439	) * &
	1440	MERGE( 1.0_wp, 0.0_wp, &
	1441	BTEST( wall_flags_0(k,j,i), 9 ) &
	1442	)
[1007]	1443	ELSE IF ( cloud_droplets ) THEN
[1342]	1444	k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
	1445	k2 = 0.61_wp * pt(k,j,i)
[2232]	1446	tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) * &
	1447	( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + &
	1448	k2 * ( q(k+1,j,i) - q(k-1,j,i) ) - &
	1449	pt(k,j,i) * ( ql(k+1,j,i) - &
	1450	ql(k-1,j,i) ) ) * dd2zu(k)&
	1451	* MERGE( 1.0_wp, 0.0_wp, &
	1452	BTEST( wall_flags_0(k,j,i), 30 ) &
	1453	) &
	1454	* MERGE( 1.0_wp, 0.0_wp, &
	1455	BTEST( wall_flags_0(k,j,i), 9 ) &
	1456	)
[1007]	1457	ENDIF
[940]	1458	ENDDO
[19]	1459
[940]	1460	IF ( use_surface_fluxes ) THEN
[2232]	1461	!
	1462	!-- Treat horizontal default surfaces, up- and downward-facing
	1463	DO l = 0, 1
	1464	surf_s = surf_def_h(l)%start_index(j,i)
	1465	surf_e = surf_def_h(l)%end_index(j,i)
	1466	DO m = surf_s, surf_e
	1467	k = surf_def_h(l)%k(m)
[1]	1468
[2232]	1469	IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
	1470	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
	1471	k2 = 0.61_wp * pt(k,j,i)
	1472	ELSE IF ( cloud_physics ) THEN
	1473	IF ( ql(k,j,i) == 0.0_wp ) THEN
	1474	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
	1475	k2 = 0.61_wp * pt(k,j,i)
	1476	ELSE
	1477	theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
	1478	temp = theta * t_d_pt(k)
	1479	k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * &
	1480	( q(k,j,i) - ql(k,j,i) ) * &
	1481	( 1.0_wp + 0.622_wp * l_d_r / temp ) ) / &
	1482	( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * &
	1483	( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
	1484	k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
	1485	ENDIF
	1486	ELSE IF ( cloud_droplets ) THEN
	1487	k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
	1488	k2 = 0.61_wp * pt(k,j,i)
	1489	ENDIF
	1490
	1491	tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * &
	1492	( k1 * surf_def_h(l)%shf(m) + &
[2329]	1493	k2 * surf_def_h(l)%qsws(m) &
	1494	) * drho_air_zw(k-1)
[2232]	1495	ENDDO
	1496	ENDDO
	1497	!
	1498	!-- Treat horizontal natural surfaces
	1499	surf_s = surf_lsm_h%start_index(j,i)
	1500	surf_e = surf_lsm_h%end_index(j,i)
	1501	DO m = surf_s, surf_e
	1502	k = surf_lsm_h%k(m)
	1503
	1504	IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
	1505	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
	1506	k2 = 0.61_wp * pt(k,j,i)
	1507	ELSE IF ( cloud_physics ) THEN
	1508	IF ( ql(k,j,i) == 0.0_wp ) THEN
	1509	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
	1510	k2 = 0.61_wp * pt(k,j,i)
	1511	ELSE
	1512	theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
	1513	temp = theta * t_d_pt(k)
	1514	k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * &
	1515	( q(k,j,i) - ql(k,j,i) ) * &
	1516	( 1.0_wp + 0.622_wp * l_d_r / temp ) ) / &
	1517	( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * &
	1518	( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
	1519	k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
	1520	ENDIF
	1521	ELSE IF ( cloud_droplets ) THEN
	1522	k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
[1342]	1523	k2 = 0.61_wp * pt(k,j,i)
[940]	1524	ENDIF
	1525
[2232]	1526	tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * &
	1527	( k1 * surf_lsm_h%shf(m) + &
[2329]	1528	k2 * surf_lsm_h%qsws(m) &
	1529	) * drho_air_zw(k-1)
[2232]	1530	ENDDO
	1531	!
	1532	!-- Treat horizontal urban surfaces
	1533	surf_s = surf_usm_h%start_index(j,i)
	1534	surf_e = surf_usm_h%end_index(j,i)
	1535	DO m = surf_s, surf_e
	1536	k = surf_usm_h%k(m)
	1537
	1538	IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
	1539	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
	1540	k2 = 0.61_wp * pt(k,j,i)
	1541	ELSE IF ( cloud_physics ) THEN
	1542	IF ( ql(k,j,i) == 0.0_wp ) THEN
	1543	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
	1544	k2 = 0.61_wp * pt(k,j,i)
	1545	ELSE
	1546	theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
	1547	temp = theta * t_d_pt(k)
	1548	k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * &
	1549	( q(k,j,i) - ql(k,j,i) ) * &
	1550	( 1.0_wp + 0.622_wp * l_d_r / temp ) ) / &
	1551	( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * &
	1552	( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
	1553	k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
	1554	ENDIF
	1555	ELSE IF ( cloud_droplets ) THEN
	1556	k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
	1557	k2 = 0.61_wp * pt(k,j,i)
	1558	ENDIF
	1559
	1560	tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * &
	1561	( k1 * surf_usm_h%shf(m) + &
[2329]	1562	k2 * surf_usm_h%qsws(m) &
	1563	) * drho_air_zw(k-1)
[2232]	1564	ENDDO
	1565
[1]	1566	ENDIF
	1567
[940]	1568	IF ( use_top_fluxes ) THEN
[2232]	1569	surf_s = surf_def_h(2)%start_index(j,i)
	1570	surf_e = surf_def_h(2)%end_index(j,i)
	1571	DO m = surf_s, surf_e
	1572	k = surf_def_h(2)%k(m)
[1]	1573
[2232]	1574
	1575
	1576	IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
[1342]	1577	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
	1578	k2 = 0.61_wp * pt(k,j,i)
[2232]	1579	ELSE IF ( cloud_physics ) THEN
	1580	IF ( ql(k,j,i) == 0.0_wp ) THEN
	1581	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
	1582	k2 = 0.61_wp * pt(k,j,i)
	1583	ELSE
	1584	theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
	1585	temp = theta * t_d_pt(k)
	1586	k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * &
	1587	( q(k,j,i) - ql(k,j,i) ) * &
	1588	( 1.0_wp + 0.622_wp * l_d_r / temp ) ) / &
	1589	( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * &
[940]	1590	( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
[2232]	1591	k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
	1592	ENDIF
	1593	ELSE IF ( cloud_droplets ) THEN
	1594	k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
	1595	k2 = 0.61_wp * pt(k,j,i)
[940]	1596	ENDIF
	1597
[2232]	1598	tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * &
	1599	( k1* surf_def_h(2)%shf(m) + &
[2329]	1600	k2 * surf_def_h(2)%qsws(m) &
	1601	) * drho_air_zw(k-1)
[2232]	1602	ENDDO
	1603
[19]	1604	ENDIF
	1605
	1606	ENDIF
	1607
[1]	1608	ENDIF
	1609
	1610	END SUBROUTINE production_e_ij
	1611
	1612
[1682]	1613	!------------------------------------------------------------------------------!
	1614	! Description:
	1615	! ------------
	1616	!> @todo Missing subroutine description.
	1617	!------------------------------------------------------------------------------!
[1]	1618	SUBROUTINE production_e_init
	1619
[1320]	1620	USE arrays_3d, &
[2329]	1621	ONLY: kh, km, drho_air_zw, u, v, zu
[1]	1622
[1320]	1623	USE control_parameters, &
[1691]	1624	ONLY: constant_flux_layer, kappa
[1320]	1625
	1626	USE indices, &
	1627	ONLY: nxl, nxlg, nxr, nxrg, nys, nysg, nyn, nyng, nzb_u_inner, &
	1628	nzb_v_inner
	1629
[2232]	1630	USE surface_mod, &
	1631	ONLY : surf_def_h, surf_def_v, surf_lsm_h, surf_usm_h
	1632
[1]	1633	IMPLICIT NONE
	1634
[2232]	1635	INTEGER(iwp) :: i !< grid index x-direction
	1636	INTEGER(iwp) :: j !< grid index y-direction
	1637	INTEGER(iwp) :: k !< grid index z-direction
	1638	INTEGER(iwp) :: l !< running index surface type (up- or downward-facing)
	1639	INTEGER(iwp) :: m !< running index surface elements
[1]	1640
[1691]	1641	IF ( constant_flux_layer ) THEN
[1]	1642	!
	1643	!-- Calculate a virtual velocity at the surface in a way that the
	1644	!-- vertical velocity gradient at k = 1 (u(k+1)-u_0) matches the
	1645	!-- Prandtl law (-w'u'/km). This gradient is used in the TKE shear
	1646	!-- production term at k=1 (see production_e_ij).
	1647	!-- The velocity gradient has to be limited in case of too small km
	1648	!-- (otherwise the timestep may be significantly reduced by large
	1649	!-- surface winds).
[106]	1650	!-- not available in case of non-cyclic boundary conditions.
[1]	1651	!-- WARNING: the exact analytical solution would require the determination
	1652	!-- of the eddy diffusivity by km = u* * kappa * zp / phi_m.
[2232]	1653	!-- Default surfaces, upward-facing
	1654	!$OMP PARALLEL DO PRIVATE(i,j,k,m)
	1655	DO m = 1, surf_def_h(0)%ns
[1]	1656
[2232]	1657	i = surf_def_h(0)%i(m)
	1658	j = surf_def_h(0)%j(m)
	1659	k = surf_def_h(0)%k(m)
	1660	!
	1661	!-- Note, calculatione of u_0 and v_0 is not fully accurate, as u/v
	1662	!-- and km are not on the same grid. Actually, a further
	1663	!-- interpolation of km onto the u/v-grid is necessary. However, the
	1664	!-- effect of this error is negligible.
	1665	surf_def_h(0)%u_0(m) = u(k+1,j,i) + surf_def_h(0)%usws(m) * &
[2329]	1666	drho_air_zw(k-1) * &
[2232]	1667	( zu(k+1) - zu(k-1) ) / &
	1668	( km(k,j,i) + 1.0E-20_wp )
	1669	surf_def_h(0)%v_0(m) = v(k+1,j,i) + surf_def_h(0)%vsws(m) * &
[2329]	1670	drho_air_zw(k-1) * &
[2232]	1671	( zu(k+1) - zu(k-1) ) / &
	1672	( km(k,j,i) + 1.0E-20_wp )
[1]	1673
[2232]	1674	IF ( ABS( u(k+1,j,i) - surf_def_h(0)%u_0(m) ) > &
	1675	ABS( u(k+1,j,i) - u(k-1,j,i) ) &
	1676	) surf_def_h(0)%u_0(m) = u(k-1,j,i)
[1]	1677
[2232]	1678	IF ( ABS( v(k+1,j,i) - surf_def_h(0)%v_0(m) ) > &
	1679	ABS( v(k+1,j,i) - v(k-1,j,i) ) &
	1680	) surf_def_h(0)%v_0(m) = v(k-1,j,i)
[1]	1681
	1682	ENDDO
[2232]	1683	!
	1684	!-- Default surfaces, downward-facing
	1685	!$OMP PARALLEL DO PRIVATE(i,j,k,m)
	1686	DO m = 1, surf_def_h(1)%ns
[1]	1687
[2232]	1688	i = surf_def_h(1)%i(m)
	1689	j = surf_def_h(1)%j(m)
	1690	k = surf_def_h(1)%k(m)
	1691	!
	1692	!-- Note, calculatione of u_0 and v_0 is not fully accurate, as u/v
	1693	!-- and km are not on the same grid. Actually, a further
	1694	!-- interpolation of km onto the u/v-grid is necessary. However, the
	1695	!-- effect of this error is negligible.
	1696	!-- In case of downward-facing surfaces, gradient is calculated
	1697	!-- between u_0 and u(k-1).
	1698	surf_def_h(1)%u_0(m) = u(k-1,j,i) - surf_def_h(1)%usws(m) * &
[2329]	1699	drho_air_zw(k-1) * &
[2232]	1700	( zu(k+1) - zu(k-1) ) / &
	1701	( km(k,j,i) + 1.0E-20_wp )
	1702	surf_def_h(1)%v_0(m) = v(k-1,j,i) - surf_def_h(1)%vsws(m) * &
[2329]	1703	drho_air_zw(k-1) * &
[2232]	1704	( zu(k+1) - zu(k-1) ) / &
	1705	( km(k,j,i) + 1.0E-20_wp )
[1]	1706
[2232]	1707	IF ( ABS( surf_def_h(1)%u_0(m) - u(k-1,j,i) ) > &
	1708	ABS( u(k+1,j,i) - u(k-1,j,i) ) &
	1709	) surf_def_h(1)%u_0(m) = u(k+1,j,i)
	1710
	1711	IF ( ABS( surf_def_h(1)%v_0(m) - v(k-1,j,i) ) > &
	1712	ABS( v(k+1,j,i) - v(k-1,j,i) ) &
	1713	) surf_def_h(1)%v_0(m) = v(k+1,j,i)
	1714
	1715	ENDDO
	1716	!
	1717	!-- Natural surfaces, upward-facing
	1718	!$OMP PARALLEL DO PRIVATE(i,j,k,m)
	1719	DO m = 1, surf_lsm_h%ns
	1720
	1721	i = surf_lsm_h%i(m)
	1722	j = surf_lsm_h%j(m)
	1723	k = surf_lsm_h%k(m)
	1724	!
	1725	!-- Note, calculatione of u_0 and v_0 is not fully accurate, as u/v
	1726	!-- and km are not on the same grid. Actually, a further
	1727	!-- interpolation of km onto the u/v-grid is necessary. However, the
	1728	!-- effect of this error is negligible.
	1729	surf_lsm_h%u_0(m) = u(k+1,j,i) + surf_lsm_h%usws(m) * &
[2329]	1730	drho_air_zw(k-1) * &
[2232]	1731	( zu(k+1) - zu(k-1) ) / &
	1732	( km(k,j,i) + 1.0E-20_wp )
	1733	surf_lsm_h%v_0(m) = v(k+1,j,i) + surf_lsm_h%vsws(m) * &
[2329]	1734	drho_air_zw(k-1) * &
[2232]	1735	( zu(k+1) - zu(k-1) ) / &
	1736	( km(k,j,i) + 1.0E-20_wp )
	1737
	1738	IF ( ABS( u(k+1,j,i) - surf_lsm_h%u_0(m) ) > &
	1739	ABS( u(k+1,j,i) - u(k-1,j,i) ) &
	1740	) surf_lsm_h%u_0(m) = u(k-1,j,i)
	1741
	1742	IF ( ABS( v(k+1,j,i) - surf_lsm_h%v_0(m) ) > &
	1743	ABS( v(k+1,j,i) - v(k-1,j,i) ) &
	1744	) surf_lsm_h%v_0(m) = v(k-1,j,i)
	1745
	1746	ENDDO
	1747	!
	1748	!-- Urban surfaces, upward-facing
	1749	!$OMP PARALLEL DO PRIVATE(i,j,k,m)
	1750	DO m = 1, surf_usm_h%ns
	1751
	1752	i = surf_usm_h%i(m)
	1753	j = surf_usm_h%j(m)
	1754	k = surf_usm_h%k(m)
	1755	!
	1756	!-- Note, calculatione of u_0 and v_0 is not fully accurate, as u/v
	1757	!-- and km are not on the same grid. Actually, a further
	1758	!-- interpolation of km onto the u/v-grid is necessary. However, the
	1759	!-- effect of this error is negligible.
	1760	surf_usm_h%u_0(m) = u(k+1,j,i) + surf_usm_h%usws(m) * &
[2329]	1761	drho_air_zw(k-1) * &
[2232]	1762	( zu(k+1) - zu(k-1) ) / &
	1763	( km(k,j,i) + 1.0E-20_wp )
	1764	surf_usm_h%v_0(m) = v(k+1,j,i) + surf_usm_h%vsws(m) * &
[2329]	1765	drho_air_zw(k-1) * &
[2232]	1766	( zu(k+1) - zu(k-1) ) / &
	1767	( km(k,j,i) + 1.0E-20_wp )
	1768
	1769	IF ( ABS( u(k+1,j,i) - surf_usm_h%u_0(m) ) > &
	1770	ABS( u(k+1,j,i) - u(k-1,j,i) ) &
	1771	) surf_usm_h%u_0(m) = u(k-1,j,i)
	1772
	1773	IF ( ABS( v(k+1,j,i) - surf_usm_h%v_0(m) ) > &
	1774	ABS( v(k+1,j,i) - v(k-1,j,i) ) &
	1775	) surf_usm_h%v_0(m) = v(k-1,j,i)
	1776
	1777	ENDDO
	1778
[1]	1779	ENDIF
	1780
	1781	END SUBROUTINE production_e_init
	1782
	1783	END MODULE production_e_mod

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |