Home

Context Navigation

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

Last change on this file since 974 was 941, checked in by raasch, 12 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 49.1 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |