Home

Context Navigation

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

Last change on this file since 222 was 208, checked in by raasch, 16 years ago
bugfix in production_e
Property svn:keywords set to `Id`
File size: 45.6 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |