Home

Context Navigation

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

Last change on this file since 418 was 410, checked in by letzel, 15 years ago
reintegrate branch letzel/masked_output into trunk; new funtionality: masked data output (not yet documented)
Property svn:keywords set to `Id`
File size: 47.4 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |