Home

Context Navigation

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

Last change on this file since 671 was 668, checked in by suehring, 14 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 47.4 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |