Home

Context Navigation

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

Last change on this file since 547 was 484, checked in by raasch, 15 years ago
typo in file headers removed
Property svn:keywords set to `Id`
File size: 47.3 KB

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