Home

Context Navigation

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

Last change on this file since 130 was 124, checked in by raasch, 17 years ago
further bugfixes for the ocean part
Property svn:keywords set to `Id`
File size: 40.4 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |