Home

Context Navigation

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

Last change on this file since 100 was 98, checked in by raasch, 17 years ago
updating comments and rc-file
Property svn:keywords set to `Id`
File size: 39.9 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |