Home

Context Navigation

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

Last change on this file since 76 was 75, checked in by raasch, 18 years ago
preliminary update for changes concerning non-cyclic boundary conditions
Property svn:keywords set to `Id`
File size: 37.4 KB

Rev	Line
[1]	1	MODULE production_e_mod
	2
	3	!------------------------------------------------------------------------------!
	4	! Actual revisions:
	5	! -----------------
[57]	6	! Wall functions now include diabatic conditions, call of routine wall_fluxes_e,
[75]	7	! reference temperature pt_reference can be used in buoyancy term,
	8	! moisture renamed humidity
[39]	9	!
	10	! Former revisions:
	11	! -----------------
	12	! $Id: production_e.f90 75 2007-03-22 09:54:05Z raasch $
	13	!
	14	! 37 2007-03-01 08:33:54Z raasch
[19]	15	! Calculation extended for gridpoint nzt, extended for given temperature /
[37]	16	! humidity fluxes at the top, wall-part is now executed in case that a
	17	! Prandtl-layer is switched on (instead of surfaces fluxes switched on)
[1]	18	!
[3]	19	! RCS Log replace by Id keyword, revision history cleaned up
	20	!
[1]	21	! Revision 1.21 2006/04/26 12:45:35 raasch
	22	! OpenMP parallelization of production_e_init
	23	!
	24	! Revision 1.1 1997/09/19 07:45:35 raasch
	25	! Initial revision
	26	!
	27	!
	28	! Description:
	29	! ------------
	30	! Production terms (shear + buoyancy) of the TKE
[37]	31	! WARNING: the case with prandtl_layer = F and use_surface_fluxes = T is
	32	! not considered well!
[1]	33	!------------------------------------------------------------------------------!
	34
[56]	35	USE wall_fluxes_mod
	36
[1]	37	PRIVATE
	38	PUBLIC production_e, production_e_init
[56]	39
[1]	40	LOGICAL, SAVE :: first_call = .TRUE.
	41
	42	REAL, DIMENSION(:,:), ALLOCATABLE, SAVE :: u_0, v_0
	43
	44	INTERFACE production_e
	45	MODULE PROCEDURE production_e
	46	MODULE PROCEDURE production_e_ij
	47	END INTERFACE production_e
	48
	49	INTERFACE production_e_init
	50	MODULE PROCEDURE production_e_init
	51	END INTERFACE production_e_init
	52
	53	CONTAINS
	54
	55
	56	!------------------------------------------------------------------------------!
	57	! Call for all grid points
	58	!------------------------------------------------------------------------------!
	59	SUBROUTINE production_e
	60
	61	USE arrays_3d
	62	USE cloud_parameters
	63	USE control_parameters
	64	USE grid_variables
	65	USE indices
	66	USE statistics
	67
	68	IMPLICIT NONE
	69
	70	INTEGER :: i, j, k
	71
	72	REAL :: def, dudx, dudy, dudz, dvdx, dvdy, dvdz, dwdx, dwdy, dwdz, &
[53]	73	k1, k2, theta, temp
[1]	74
[56]	75	! REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: usvs, vsus, wsus, wsvs
	76	REAL, DIMENSION(nzb:nzt+1) :: usvs, vsus, wsus, wsvs
[1]	77
[56]	78	!
	79	!-- First calculate horizontal momentum flux u'v', w'v', v'u', w'u' at
	80	!-- vertical walls, if neccessary
	81	!-- So far, results are slightly different from the ij-Version.
	82	!-- Therefore, ij-Version is called further below within the ij-loops.
	83	! IF ( topography /= 'flat' ) THEN
	84	! CALL wall_fluxes_e( usvs, 1.0, 0.0, 0.0, 0.0, wall_e_y )
	85	! CALL wall_fluxes_e( wsvs, 0.0, 0.0, 1.0, 0.0, wall_e_y )
	86	! CALL wall_fluxes_e( vsus, 0.0, 1.0, 0.0, 0.0, wall_e_x )
	87	! CALL wall_fluxes_e( wsus, 0.0, 0.0, 0.0, 1.0, wall_e_x )
	88	! ENDIF
[53]	89
[1]	90	!
	91	!-- Calculate TKE production by shear
	92	DO i = nxl, nxr
	93
	94	DO j = nys, nyn
[19]	95	DO k = nzb_diff_s_outer(j,i), nzt
[1]	96
	97	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
	98	dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
	99	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
	100	dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
	101	u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
	102
	103	dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
	104	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
	105	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
	106	dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
	107	v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
	108
	109	dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
	110	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
	111	dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
	112	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
	113	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
	114
	115	def = 2.0 * ( dudx2 + dvdy2 + dwdz**2 ) + &
	116	dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz**2 + &
	117	dvdz*2 + 2.0 ( dvdxdudy + dwdxdudz + dwdy*dvdz )
	118
	119	IF ( def < 0.0 ) def = 0.0
	120
	121	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
	122
	123	ENDDO
	124	ENDDO
	125
[37]	126	IF ( prandtl_layer ) THEN
[1]	127
	128	!
[55]	129	!-- Position beneath wall
	130	!-- (2) - Will allways be executed.
	131	!-- 'bottom and wall: use u_0,v_0 and wall functions'
[1]	132	DO j = nys, nyn
	133
	134	IF ( ( wall_e_x(j,i) /= 0.0 ) .OR. ( wall_e_y(j,i) /= 0.0 ) ) &
	135	THEN
	136
	137	k = nzb_diff_s_inner(j,i) - 1
	138	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
[53]	139	dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
	140	u_0(j,i) - u_0(j,i+1) ) * dd2zu(k)
	141	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
	142	dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
	143	v_0(j,i) - v_0(j+1,i) ) * dd2zu(k)
	144	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
	145
[1]	146	IF ( wall_e_y(j,i) /= 0.0 ) THEN
[53]	147	CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
	148	usvs, 1.0, 0.0, 0.0, 0.0 )
	149	dudy = wall_e_y(j,i) * usvs(k) / km(k,j,i)
[56]	150	! dudy = wall_e_y(j,i) * usvs(k,j,i) / km(k,j,i)
[53]	151	CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
	152	wsvs, 0.0, 0.0, 1.0, 0.0 )
	153	dwdy = wall_e_y(j,i) * wsvs(k) / km(k,j,i)
[56]	154	! dwdy = wall_e_y(j,i) * wsvs(k,j,i) / km(k,j,i)
[1]	155	ELSE
	156	dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
	157	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
[53]	158	dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
	159	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
[1]	160	ENDIF
	161
	162	IF ( wall_e_x(j,i) /= 0.0 ) THEN
[53]	163	CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
	164	vsus, 0.0, 1.0, 0.0, 0.0 )
	165	dvdx = wall_e_x(j,i) * vsus(k) / km(k,j,i)
[56]	166	! dvdx = wall_e_x(j,i) * vsus(k,j,i) / km(k,j,i)
[53]	167	CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
	168	wsus, 0.0, 0.0, 0.0, 1.0 )
	169	dwdx = wall_e_x(j,i) * wsus(k) / km(k,j,i)
[56]	170	! dwdx = wall_e_x(j,i) * wsus(k,j,i) / km(k,j,i)
[1]	171	ELSE
	172	dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
	173	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
	174	dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
	175	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
	176	ENDIF
	177
	178	def = 2.0 * ( dudx2 + dvdy2 + dwdz**2 ) + &
	179	dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz**2 + &
	180	dvdz*2 + 2.0 ( dvdxdudy + dwdxdudz + dwdy*dvdz )
	181
	182	IF ( def < 0.0 ) def = 0.0
	183
	184	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
	185
	186
	187	!
[55]	188	!-- (3) - will be executed only, if there is at least one level
	189	!-- between (2) and (4), i.e. the topography must have a
	190	!-- minimum height of 2 dz. Wall fluxes for this case have
	191	!-- already been calculated for (2).
	192	!-- 'wall only: use wall functions'
[1]	193
	194	DO k = nzb_diff_s_inner(j,i), nzb_diff_s_outer(j,i)-2
	195
	196	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
[53]	197	dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
	198	u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
	199	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
	200	dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
	201	v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
	202	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
	203
[1]	204	IF ( wall_e_y(j,i) /= 0.0 ) THEN
[53]	205	dudy = wall_e_y(j,i) * usvs(k) / km(k,j,i)
[56]	206	! dudy = wall_e_y(j,i) * usvs(k,j,i) / km(k,j,i)
[53]	207	dwdy = wall_e_y(j,i) * wsvs(k) / km(k,j,i)
[56]	208	! dwdy = wall_e_y(j,i) * wsvs(k,j,i) / km(k,j,i)
[1]	209	ELSE
	210	dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
	211	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
[53]	212	dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
	213	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
[1]	214	ENDIF
	215
	216	IF ( wall_e_x(j,i) /= 0.0 ) THEN
[53]	217	dvdx = wall_e_x(j,i) * vsus(k) / km(k,j,i)
[56]	218	! dvdx = wall_e_x(j,i) * vsus(k,j,i) / km(k,j,i)
[53]	219	dwdx = wall_e_x(j,i) * wsus(k) / km(k,j,i)
[56]	220	! dwdx = wall_e_x(j,i) * wsus(k,j,i) / km(k,j,i)
[1]	221	ELSE
	222	dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
	223	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
	224	dwdx = 0.25 * ( 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	ENDIF
	227
	228	def = 2.0 * ( dudx2 + dvdy2 + dwdz**2 ) + &
	229	dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz**2 + &
	230	dvdz*2 + 2.0 ( dvdxdudy + dwdxdudz + dwdy*dvdz )
	231
	232	IF ( def < 0.0 ) def = 0.0
	233
	234	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
	235
	236	ENDDO
	237
	238	ENDIF
	239
	240	ENDDO
	241
	242	!
[55]	243	!-- (4) - will allways be executed.
	244	!-- 'special case: free atmosphere' (as for case (0))
[1]	245	DO j = nys, nyn
	246
	247	IF ( ( wall_e_x(j,i) /= 0.0 ) .OR. ( wall_e_y(j,i) /= 0.0 ) ) &
	248	THEN
	249
	250	k = nzb_diff_s_outer(j,i)-1
	251
	252	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
	253	dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
	254	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
	255	dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
	256	u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
	257
	258	dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
	259	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
	260	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
	261	dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
	262	v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
	263
	264	dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
	265	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
	266	dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
	267	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
	268	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
	269
	270	def = 2.0 * ( dudx2 + dvdy2 + dwdz**2 ) + &
	271	dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz**2 + &
	272	dvdz*2 + 2.0 ( dvdxdudy + dwdxdudz + dwdy*dvdz )
	273
	274	IF ( def < 0.0 ) def = 0.0
	275
	276	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
	277
	278	ENDIF
	279
	280	ENDDO
	281
	282	!
[55]	283	!-- Position without adjacent wall
	284	!-- (1) - will allways be executed.
	285	!-- 'bottom only: use u_0,v_0'
[1]	286	DO j = nys, nyn
	287
	288	IF ( ( wall_e_x(j,i) == 0.0 ) .AND. ( wall_e_y(j,i) == 0.0 ) ) &
	289	THEN
	290
	291	k = nzb_diff_s_inner(j,i)-1
	292
	293	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
	294	dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
	295	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
	296	dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
	297	u_0(j,i) - u_0(j,i+1) ) * dd2zu(k)
	298
	299	dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
	300	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
	301	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
	302	dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
	303	v_0(j,i) - v_0(j+1,i) ) * dd2zu(k)
	304
	305	dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
	306	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
	307	dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
	308	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
	309	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
	310
	311	def = 2.0 * ( dudx2 + dvdy2 + dwdz**2 ) + &
	312	dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz**2 + &
	313	dvdz*2 + 2.0 ( dvdxdudy + dwdxdudz + dwdy*dvdz )
	314
	315	IF ( def < 0.0 ) def = 0.0
	316
	317	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
	318
	319	ENDIF
	320
	321	ENDDO
	322
[37]	323	ELSEIF ( use_surface_fluxes ) THEN
	324
	325	DO j = nys, nyn
	326
	327	k = nzb_diff_s_outer(j,i)-1
	328
	329	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
	330	dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
	331	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
	332	dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
	333	u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
	334
	335	dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
	336	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
	337	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
	338	dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
	339	v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
	340
	341	dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
	342	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
	343	dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
	344	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
	345	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
	346
	347	def = 2.0 * ( dudx2 + dvdy2 + dwdz**2 ) + &
	348	dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz**2 + &
	349	dvdz*2 + 2.0 ( dvdxdudy + dwdxdudz + dwdy*dvdz )
	350
	351	IF ( def < 0.0 ) def = 0.0
	352
	353	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
	354
	355	ENDDO
	356
[1]	357	ENDIF
	358
	359	!
	360	!-- Calculate TKE production by buoyancy
[75]	361	IF ( .NOT. humidity ) THEN
[1]	362
[57]	363	IF ( use_pt_reference ) THEN
[1]	364
[57]	365	DO j = nys, nyn
	366	DO k = nzb_diff_s_inner(j,i), nzt_diff
	367	tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g/pt_reference * &
	368	( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k)
	369	ENDDO
	370
	371	IF ( use_surface_fluxes ) THEN
	372	k = nzb_diff_s_inner(j,i)-1
	373	tend(k,j,i) = tend(k,j,i) + g / pt_reference * shf(j,i)
	374	ENDIF
	375
	376	IF ( use_top_fluxes ) THEN
	377	k = nzt
	378	tend(k,j,i) = tend(k,j,i) + g / pt_reference * tswst(j,i)
	379	ENDIF
[1]	380	ENDDO
[19]	381
[57]	382	ELSE
[1]	383
[57]	384	DO j = nys, nyn
	385	DO k = nzb_diff_s_inner(j,i), nzt_diff
	386	tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / pt(k,j,i) * &
	387	( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k)
	388	ENDDO
[19]	389
[57]	390	IF ( use_surface_fluxes ) THEN
	391	k = nzb_diff_s_inner(j,i)-1
	392	tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * shf(j,i)
	393	ENDIF
[1]	394
[57]	395	IF ( use_top_fluxes ) THEN
	396	k = nzt
	397	tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * tswst(j,i)
	398	ENDIF
	399	ENDDO
	400
	401	ENDIF
	402
[1]	403	ELSE
	404
	405	DO j = nys, nyn
	406
[19]	407	DO k = nzb_diff_s_inner(j,i), nzt_diff
[1]	408
	409	IF ( .NOT. cloud_physics ) THEN
	410	k1 = 1.0 + 0.61 * q(k,j,i)
	411	k2 = 0.61 * pt(k,j,i)
	412	ELSE
	413	IF ( ql(k,j,i) == 0.0 ) THEN
	414	k1 = 1.0 + 0.61 * q(k,j,i)
	415	k2 = 0.61 * pt(k,j,i)
	416	ELSE
	417	theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
	418	temp = theta * t_d_pt(k)
	419	k1 = ( 1.0 - q(k,j,i) + 1.61 * &
	420	( q(k,j,i) - ql(k,j,i) ) * &
	421	( 1.0 + 0.622 * l_d_r / temp ) ) / &
	422	( 1.0 + 0.622 * l_d_r * l_d_cp * &
	423	( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
	424	k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
	425	ENDIF
	426	ENDIF
	427
	428	tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) * &
	429	( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + &
	430	k2 * ( q(k+1,j,i) - q(k-1,j,i) ) &
	431	) * dd2zu(k)
	432	ENDDO
	433
	434	ENDDO
	435
	436	IF ( use_surface_fluxes ) THEN
	437
	438	DO j = nys, nyn
	439
[53]	440	k = nzb_diff_s_inner(j,i)-1
[1]	441
	442	IF ( .NOT. cloud_physics ) THEN
	443	k1 = 1.0 + 0.61 * q(k,j,i)
	444	k2 = 0.61 * pt(k,j,i)
	445	ELSE
	446	IF ( ql(k,j,i) == 0.0 ) THEN
	447	k1 = 1.0 + 0.61 * q(k,j,i)
	448	k2 = 0.61 * pt(k,j,i)
	449	ELSE
	450	theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
	451	temp = theta * t_d_pt(k)
	452	k1 = ( 1.0 - q(k,j,i) + 1.61 * &
	453	( q(k,j,i) - ql(k,j,i) ) * &
	454	( 1.0 + 0.622 * l_d_r / temp ) ) / &
	455	( 1.0 + 0.622 * l_d_r * l_d_cp * &
	456	( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
	457	k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
	458	ENDIF
	459	ENDIF
	460
	461	tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * &
	462	( k1* shf(j,i) + k2 * qsws(j,i) )
	463	ENDDO
	464
	465	ENDIF
	466
[19]	467	IF ( use_top_fluxes ) THEN
	468
	469	DO j = nys, nyn
	470
	471	k = nzt
	472
	473	IF ( .NOT. cloud_physics ) THEN
	474	k1 = 1.0 + 0.61 * q(k,j,i)
	475	k2 = 0.61 * pt(k,j,i)
	476	ELSE
	477	IF ( ql(k,j,i) == 0.0 ) THEN
	478	k1 = 1.0 + 0.61 * q(k,j,i)
	479	k2 = 0.61 * pt(k,j,i)
	480	ELSE
	481	theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
	482	temp = theta * t_d_pt(k)
	483	k1 = ( 1.0 - q(k,j,i) + 1.61 * &
	484	( q(k,j,i) - ql(k,j,i) ) * &
	485	( 1.0 + 0.622 * l_d_r / temp ) ) / &
	486	( 1.0 + 0.622 * l_d_r * l_d_cp * &
	487	( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
	488	k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
	489	ENDIF
	490	ENDIF
	491
	492	tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * &
	493	( k1* tswst(j,i) + k2 * qswst(j,i) )
	494	ENDDO
	495
	496	ENDIF
	497
[1]	498	ENDIF
	499
	500	ENDDO
	501
	502	END SUBROUTINE production_e
	503
	504
	505	!------------------------------------------------------------------------------!
	506	! Call for grid point i,j
	507	!------------------------------------------------------------------------------!
	508	SUBROUTINE production_e_ij( i, j )
	509
	510	USE arrays_3d
	511	USE cloud_parameters
	512	USE control_parameters
	513	USE grid_variables
	514	USE indices
	515	USE statistics
	516
	517	IMPLICIT NONE
	518
	519	INTEGER :: i, j, k
	520
	521	REAL :: def, dudx, dudy, dudz, dvdx, dvdy, dvdz, dwdx, dwdy, dwdz, &
[53]	522	k1, k2, theta, temp
[1]	523
[56]	524	REAL, DIMENSION(nzb:nzt+1) :: usvs, vsus, wsus, wsvs
[53]	525
[1]	526	!
	527	!-- Calculate TKE production by shear
[19]	528	DO k = nzb_diff_s_outer(j,i), nzt
[1]	529
	530	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
	531	dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
	532	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
	533	dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
	534	u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
	535
	536	dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
	537	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
	538	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
	539	dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
	540	v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
	541
	542	dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
	543	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
	544	dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
	545	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
	546	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
	547
	548	def = 2.0 * ( dudx2 + dvdy2 + dwdz**2 ) &
	549	+ dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz2 + dvdz2 &
	550	+ 2.0 * ( dvdxdudy + dwdxdudz + dwdy*dvdz )
	551
	552	IF ( def < 0.0 ) def = 0.0
	553
	554	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
	555
	556	ENDDO
	557
[37]	558	IF ( prandtl_layer ) THEN
[1]	559
	560	IF ( ( wall_e_x(j,i) /= 0.0 ) .OR. ( wall_e_y(j,i) /= 0.0 ) ) THEN
[55]	561
[1]	562	!
[55]	563	!-- Position beneath wall
	564	!-- (2) - Will allways be executed.
	565	!-- 'bottom and wall: use u_0,v_0 and wall functions'
[1]	566	k = nzb_diff_s_inner(j,i)-1
	567
	568	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
[53]	569	dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
	570	u_0(j,i) - u_0(j,i+1) ) * dd2zu(k)
	571	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
	572	dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
	573	v_0(j,i) - v_0(j+1,i) ) * dd2zu(k)
	574	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
	575
[1]	576	IF ( wall_e_y(j,i) /= 0.0 ) THEN
[53]	577	CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
	578	usvs, 1.0, 0.0, 0.0, 0.0 )
	579	dudy = wall_e_y(j,i) * usvs(k) / km(k,j,i)
	580	CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
	581	wsvs, 0.0, 0.0, 1.0, 0.0 )
	582	dwdy = wall_e_y(j,i) * wsvs(k) / km(k,j,i)
[1]	583	ELSE
	584	dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
	585	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
[53]	586	dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
	587	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
[1]	588	ENDIF
	589
	590	IF ( wall_e_x(j,i) /= 0.0 ) THEN
[53]	591	CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
	592	vsus, 0.0, 1.0, 0.0, 0.0 )
	593	dvdx = wall_e_x(j,i) * vsus(k) / km(k,j,i)
	594	CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
	595	wsus, 0.0, 0.0, 0.0, 1.0 )
	596	dwdx = wall_e_x(j,i) * wsus(k) / km(k,j,i)
[1]	597	ELSE
	598	dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
	599	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
	600	dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
	601	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
	602	ENDIF
	603
	604	def = 2.0 * ( dudx2 + dvdy2 + dwdz**2 ) + &
	605	dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz**2 + &
	606	dvdz*2 + 2.0 ( dvdxdudy + dwdxdudz + dwdy*dvdz )
	607
	608	IF ( def < 0.0 ) def = 0.0
	609
	610	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
	611
	612	!
[55]	613	!-- (3) - will be executed only, if there is at least one level
	614	!-- between (2) and (4), i.e. the topography must have a
	615	!-- minimum height of 2 dz. Wall fluxes for this case have
	616	!-- already been calculated for (2).
	617	!-- 'wall only: use wall functions'
[1]	618	DO k = nzb_diff_s_inner(j,i), nzb_diff_s_outer(j,i)-2
	619
	620	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
[53]	621	dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
	622	u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
	623	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
	624	dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
	625	v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
	626	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
	627
[1]	628	IF ( wall_e_y(j,i) /= 0.0 ) THEN
[53]	629	dudy = wall_e_y(j,i) * usvs(k) / km(k,j,i)
	630	dwdy = wall_e_y(j,i) * wsvs(k) / km(k,j,i)
[1]	631	ELSE
	632	dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
	633	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
[53]	634	dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
	635	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
[1]	636	ENDIF
	637
	638	IF ( wall_e_x(j,i) /= 0.0 ) THEN
[53]	639	dvdx = wall_e_x(j,i) * vsus(k) / km(k,j,i)
	640	dwdx = wall_e_x(j,i) * wsus(k) / km(k,j,i)
[1]	641	ELSE
	642	dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
	643	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
	644	dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
	645	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
	646	ENDIF
	647
	648	def = 2.0 * ( dudx2 + dvdy2 + dwdz**2 ) + &
	649	dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz**2 + &
	650	dvdz*2 + 2.0 ( dvdxdudy + dwdxdudz + dwdy*dvdz )
	651
	652	IF ( def < 0.0 ) def = 0.0
	653
	654	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
	655
	656	ENDDO
	657
	658	!
[55]	659	!-- (4) - will allways be executed.
	660	!-- 'special case: free atmosphere' (as for case (0))
[1]	661	k = nzb_diff_s_outer(j,i)-1
	662
	663	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
	664	dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
	665	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
	666	dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
	667	u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
	668
	669	dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
	670	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
	671	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
	672	dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
	673	v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
	674
	675	dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
	676	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
	677	dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
	678	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
	679	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
	680
	681	def = 2.0 * ( dudx2 + dvdy2 + dwdz**2 ) + &
	682	dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz**2 + &
	683	dvdz*2 + 2.0 ( dvdxdudy + dwdxdudz + dwdy*dvdz )
	684
	685	IF ( def < 0.0 ) def = 0.0
	686
	687	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
	688
	689	ELSE
	690
	691	!
[55]	692	!-- Position without adjacent wall
	693	!-- (1) - will allways be executed.
	694	!-- 'bottom only: use u_0,v_0'
[1]	695	k = nzb_diff_s_inner(j,i)-1
	696
	697	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
	698	dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
	699	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
	700	dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
	701	u_0(j,i) - u_0(j,i+1) ) * dd2zu(k)
	702
	703	dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
	704	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
	705	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
	706	dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
	707	v_0(j,i) - v_0(j+1,i) ) * dd2zu(k)
	708
	709	dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
	710	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
	711	dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
	712	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
	713	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
	714
	715	def = 2.0 * ( dudx2 + dvdy2 + dwdz**2 ) &
	716	+ dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz2 + dvdz2 &
	717	+ 2.0 * ( dvdxdudy + dwdxdudz + dwdy*dvdz )
	718
	719	IF ( def < 0.0 ) def = 0.0
	720
	721	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
	722
	723	ENDIF
	724
[37]	725	ELSEIF ( use_surface_fluxes ) THEN
	726
	727	k = nzb_diff_s_outer(j,i)-1
	728
	729	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
	730	dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
	731	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
	732	dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
	733	u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
	734
	735	dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
	736	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
	737	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
	738	dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
	739	v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
	740
	741	dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
	742	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
	743	dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
	744	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
	745	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
	746
	747	def = 2.0 * ( dudx2 + dvdy2 + dwdz**2 ) + &
	748	dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz**2 + &
	749	dvdz*2 + 2.0 ( dvdxdudy + dwdxdudz + dwdy*dvdz )
	750
	751	IF ( def < 0.0 ) def = 0.0
	752
	753	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
	754
[1]	755	ENDIF
	756
	757	!
	758	!-- Calculate TKE production by buoyancy
[75]	759	IF ( .NOT. humidity ) THEN
[1]	760
[57]	761	IF ( use_pt_reference ) THEN
[19]	762
[57]	763	DO k = nzb_diff_s_inner(j,i), nzt_diff
	764	tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / pt_reference * &
	765	( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k)
	766	ENDDO
[1]	767
[57]	768	IF ( use_surface_fluxes ) THEN
	769	k = nzb_diff_s_inner(j,i)-1
	770	tend(k,j,i) = tend(k,j,i) + g / pt_reference * shf(j,i)
	771	ENDIF
[19]	772
[57]	773	IF ( use_top_fluxes ) THEN
	774	k = nzt
	775	tend(k,j,i) = tend(k,j,i) + g / pt_reference * tswst(j,i)
	776	ENDIF
	777
	778	ELSE
	779
	780	DO k = nzb_diff_s_inner(j,i), nzt_diff
	781	tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / pt(k,j,i) * &
	782	( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k)
	783	ENDDO
	784
	785	IF ( use_surface_fluxes ) THEN
	786	k = nzb_diff_s_inner(j,i)-1
	787	tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * shf(j,i)
	788	ENDIF
	789
	790	IF ( use_top_fluxes ) THEN
	791	k = nzt
	792	tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * tswst(j,i)
	793	ENDIF
	794
	795	ENDIF
	796
[1]	797	ELSE
	798
[19]	799	DO k = nzb_diff_s_inner(j,i), nzt_diff
[1]	800
	801	IF ( .NOT. cloud_physics ) THEN
	802	k1 = 1.0 + 0.61 * q(k,j,i)
	803	k2 = 0.61 * pt(k,j,i)
	804	ELSE
	805	IF ( ql(k,j,i) == 0.0 ) THEN
	806	k1 = 1.0 + 0.61 * q(k,j,i)
	807	k2 = 0.61 * pt(k,j,i)
	808	ELSE
	809	theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
	810	temp = theta * t_d_pt(k)
	811	k1 = ( 1.0 - q(k,j,i) + 1.61 * &
	812	( q(k,j,i) - ql(k,j,i) ) * &
	813	( 1.0 + 0.622 * l_d_r / temp ) ) / &
	814	( 1.0 + 0.622 * l_d_r * l_d_cp * &
	815	( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
	816	k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
	817	ENDIF
	818	ENDIF
	819
	820	tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) * &
	821	( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + &
	822	k2 * ( q(k+1,j,i) - q(k-1,j,i) ) &
	823	) * dd2zu(k)
	824	ENDDO
[19]	825
[1]	826	IF ( use_surface_fluxes ) THEN
	827	k = nzb_diff_s_inner(j,i)-1
	828
	829	IF ( .NOT. cloud_physics ) THEN
	830	k1 = 1.0 + 0.61 * q(k,j,i)
	831	k2 = 0.61 * pt(k,j,i)
	832	ELSE
	833	IF ( ql(k,j,i) == 0.0 ) THEN
	834	k1 = 1.0 + 0.61 * q(k,j,i)
	835	k2 = 0.61 * pt(k,j,i)
	836	ELSE
	837	theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
	838	temp = theta * t_d_pt(k)
	839	k1 = ( 1.0 - q(k,j,i) + 1.61 * &
	840	( q(k,j,i) - ql(k,j,i) ) * &
	841	( 1.0 + 0.622 * l_d_r / temp ) ) / &
	842	( 1.0 + 0.622 * l_d_r * l_d_cp * &
	843	( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
	844	k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
	845	ENDIF
	846	ENDIF
	847
	848	tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * &
	849	( k1* shf(j,i) + k2 * qsws(j,i) )
	850	ENDIF
	851
[19]	852	IF ( use_top_fluxes ) THEN
	853	k = nzt
	854
	855	IF ( .NOT. cloud_physics ) THEN
	856	k1 = 1.0 + 0.61 * q(k,j,i)
	857	k2 = 0.61 * pt(k,j,i)
	858	ELSE
	859	IF ( ql(k,j,i) == 0.0 ) THEN
	860	k1 = 1.0 + 0.61 * q(k,j,i)
	861	k2 = 0.61 * pt(k,j,i)
	862	ELSE
	863	theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
	864	temp = theta * t_d_pt(k)
	865	k1 = ( 1.0 - q(k,j,i) + 1.61 * &
	866	( q(k,j,i) - ql(k,j,i) ) * &
	867	( 1.0 + 0.622 * l_d_r / temp ) ) / &
	868	( 1.0 + 0.622 * l_d_r * l_d_cp * &
	869	( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
	870	k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
	871	ENDIF
	872	ENDIF
	873
	874	tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * &
	875	( k1* tswst(j,i) + k2 * qswst(j,i) )
	876	ENDIF
	877
[1]	878	ENDIF
	879
	880	END SUBROUTINE production_e_ij
	881
	882
	883	SUBROUTINE production_e_init
	884
	885	USE arrays_3d
	886	USE control_parameters
	887	USE grid_variables
	888	USE indices
	889
	890	IMPLICIT NONE
	891
	892	INTEGER :: i, j, ku, kv
	893
[37]	894	IF ( prandtl_layer ) THEN
[1]	895
	896	IF ( first_call ) THEN
	897	ALLOCATE( u_0(nys-1:nyn+1,nxl-1:nxr+1), &
	898	v_0(nys-1:nyn+1,nxl-1:nxr+1) )
	899	first_call = .FALSE.
	900	ENDIF
	901
	902	!
	903	!-- Calculate a virtual velocity at the surface in a way that the
	904	!-- vertical velocity gradient at k = 1 (u(k+1)-u_0) matches the
	905	!-- Prandtl law (-w'u'/km). This gradient is used in the TKE shear
	906	!-- production term at k=1 (see production_e_ij).
	907	!-- The velocity gradient has to be limited in case of too small km
	908	!-- (otherwise the timestep may be significantly reduced by large
	909	!-- surface winds).
	910	!-- WARNING: the exact analytical solution would require the determination
	911	!-- of the eddy diffusivity by km = u* * kappa * zp / phi_m.
	912	!$OMP PARALLEL DO PRIVATE( ku, kv )
[73]	913	DO i = nxl, nxr
	914	DO j = nys, nyn
[1]	915
	916	ku = nzb_u_inner(j,i)+1
	917	kv = nzb_v_inner(j,i)+1
	918
	919	u_0(j,i) = u(ku+1,j,i) + usws(j,i) * ( zu(ku+1) - zu(ku-1) ) / &
	920	( 0.5 * ( km(ku,j,i) + km(ku,j,i-1) ) + &
	921	1.0E-20 )
	922	! ( us(j,i) * kappa * zu(1) )
	923	v_0(j,i) = v(kv+1,j,i) + vsws(j,i) * ( zu(kv+1) - zu(kv-1) ) / &
	924	( 0.5 * ( km(kv,j,i) + km(kv,j-1,i) ) + &
	925	1.0E-20 )
	926	! ( us(j,i) * kappa * zu(1) )
	927
	928	IF ( ABS( u(ku+1,j,i) - u_0(j,i) ) > &
	929	ABS( u(ku+1,j,i) - u(ku-1,j,i) ) ) u_0(j,i) = u(ku-1,j,i)
	930	IF ( ABS( v(kv+1,j,i) - v_0(j,i) ) > &
	931	ABS( v(kv+1,j,i) - v(kv-1,j,i) ) ) v_0(j,i) = v(kv-1,j,i)
	932
	933	ENDDO
	934	ENDDO
	935
	936	CALL exchange_horiz_2d( u_0 )
	937	CALL exchange_horiz_2d( v_0 )
	938
	939	ENDIF
	940
	941	END SUBROUTINE production_e_init
	942
	943	END MODULE production_e_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |