Home

Context Navigation

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

Last change on this file since 47 was 39, checked in by raasch, 18 years ago
comments prepared for 3.1c
Property svn:keywords set to `Id`
File size: 36.6 KB

Rev	Line
[1]	1	MODULE production_e_mod
	2
	3	!------------------------------------------------------------------------------!
	4	! Actual revisions:
	5	! -----------------
[39]	6	!
	7	!
	8	! Former revisions:
	9	! -----------------
	10	! $Id: production_e.f90 39 2007-03-01 12:46:59Z raasch $
	11	!
	12	! 37 2007-03-01 08:33:54Z raasch
[19]	13	! Calculation extended for gridpoint nzt, extended for given temperature /
[37]	14	! humidity fluxes at the top, wall-part is now executed in case that a
	15	! Prandtl-layer is switched on (instead of surfaces fluxes switched on)
[1]	16	!
[3]	17	! RCS Log replace by Id keyword, revision history cleaned up
	18	!
[1]	19	! Revision 1.21 2006/04/26 12:45:35 raasch
	20	! OpenMP parallelization of production_e_init
	21	!
	22	! Revision 1.1 1997/09/19 07:45:35 raasch
	23	! Initial revision
	24	!
	25	!
	26	! Description:
	27	! ------------
	28	! Production terms (shear + buoyancy) of the TKE
[37]	29	! WARNING: the case with prandtl_layer = F and use_surface_fluxes = T is
	30	! not considered well!
[1]	31	!------------------------------------------------------------------------------!
	32
	33	PRIVATE
	34	PUBLIC production_e, production_e_init
	35
	36	LOGICAL, SAVE :: first_call = .TRUE.
	37
	38	REAL, DIMENSION(:,:), ALLOCATABLE, SAVE :: u_0, v_0
	39
	40	INTERFACE production_e
	41	MODULE PROCEDURE production_e
	42	MODULE PROCEDURE production_e_ij
	43	END INTERFACE production_e
	44
	45	INTERFACE production_e_init
	46	MODULE PROCEDURE production_e_init
	47	END INTERFACE production_e_init
	48
	49	CONTAINS
	50
	51
	52	!------------------------------------------------------------------------------!
	53	! Call for all grid points
	54	!------------------------------------------------------------------------------!
	55	SUBROUTINE production_e
	56
	57	USE arrays_3d
	58	USE cloud_parameters
	59	USE control_parameters
	60	USE grid_variables
	61	USE indices
	62	USE statistics
	63
	64	IMPLICIT NONE
	65
	66	INTEGER :: i, j, k
	67
	68	REAL :: def, dudx, dudy, dudz, dvdx, dvdy, dvdz, dwdx, dwdy, dwdz, &
	69	k1, k2, theta, temp, usvs, vsus, wsus, wsvs
	70
	71
	72	!
	73	!-- Calculate TKE production by shear
	74	DO i = nxl, nxr
	75
	76	DO j = nys, nyn
[19]	77	DO k = nzb_diff_s_outer(j,i), nzt
[1]	78
	79	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
	80	dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
	81	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
	82	dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
	83	u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
	84
	85	dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
	86	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
	87	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
	88	dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
	89	v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
	90
	91	dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
	92	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
	93	dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
	94	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
	95	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
	96
	97	def = 2.0 * ( dudx2 + dvdy2 + dwdz**2 ) + &
	98	dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz**2 + &
	99	dvdz*2 + 2.0 ( dvdxdudy + dwdxdudz + dwdy*dvdz )
	100
	101	IF ( def < 0.0 ) def = 0.0
	102
	103	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
	104
	105	ENDDO
	106	ENDDO
	107
[37]	108	IF ( prandtl_layer ) THEN
[1]	109
	110	!
	111	!-- Position neben Gebaeudewand
	112	!-- 2 - Wird immer ausgefuehrt. 'Boden und Wand:
	113	!-- u_0,v_0 und Wall functions'
	114	DO j = nys, nyn
	115
	116	IF ( ( wall_e_x(j,i) /= 0.0 ) .OR. ( wall_e_y(j,i) /= 0.0 ) ) &
	117	THEN
	118
	119	k = nzb_diff_s_inner(j,i) - 1
	120	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
	121	IF ( wall_e_y(j,i) /= 0.0 ) THEN
	122	usvs = kappa * 0.5 * ( u(k,j,i) + u(k,j,i+1) ) &
	123	/ LOG( 0.5 * dy / z0(j,i) )
	124	usvs = usvs * ABS( usvs )
	125	dudy = wall_e_y(j,i) * usvs / km(k,j,i)
	126	ELSE
	127	dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
	128	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
	129	ENDIF
	130	dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
	131	u_0(j,i) - u_0(j,i+1) ) * dd2zu(k)
	132
	133	IF ( wall_e_x(j,i) /= 0.0 ) THEN
	134	vsus = kappa * 0.5 * ( v(k,j,i) + v(k,j+1,i) ) &
	135	/ LOG( 0.5 * dx / z0(j,i))
	136	vsus = vsus * ABS( vsus )
	137	dvdx = wall_e_x(j,i) * vsus / km(k,j,i)
	138	ELSE
	139	dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
	140	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
	141	ENDIF
	142	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
	143	dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
	144	v_0(j,i) - v_0(j+1,i) ) * dd2zu(k)
	145
	146	IF ( wall_e_x(j,i) /= 0.0 ) THEN
	147	wsus = kappa * 0.5 * ( w(k,j,i) + w(k-1,j,i) ) &
	148	/ LOG( 0.5 * dx / z0(j,i))
	149	wsus = wsus * ABS( wsus )
	150	dwdx = wall_e_x(j,i) * wsus / km(k,j,i)
	151	ELSE
	152	dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
	153	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
	154	ENDIF
	155	IF ( wall_e_y(j,i) /= 0.0 ) THEN
	156	wsvs = kappa * ( w(k,j,i) + w(k-1,j,i) ) &
	157	/ LOG( 0.5 * dy / z0(j,i))
	158	wsvs = wsvs * ABS( wsvs )
	159	dwdy = wall_e_y(j,i) * wsvs / km(k,j,i)
	160	ELSE
	161	dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
	162	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
	163	ENDIF
	164	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
	165
	166	def = 2.0 * ( dudx2 + dvdy2 + dwdz**2 ) + &
	167	dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz**2 + &
	168	dvdz*2 + 2.0 ( dvdxdudy + dwdxdudz + dwdy*dvdz )
	169
	170	IF ( def < 0.0 ) def = 0.0
	171
	172	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
	173
	174	ENDIF
	175
	176	ENDDO
	177
	178	!
	179	!-- 3 - Wird nur ausgefuehrt, wenn mindestens ein Niveau
	180	!-- zwischen 2 und 4 liegt, d.h. ab einer Gebaeudemindest-
	181	!-- hoehe von 2 dz. 'Nur Wand: Wall functions'
	182	DO j = nys, nyn
	183
	184	IF ( ( wall_e_x(j,i) /= 0.0 ) .OR. ( wall_e_y(j,i) /= 0.0 ) ) &
	185	THEN
	186
	187	DO k = nzb_diff_s_inner(j,i), nzb_diff_s_outer(j,i)-2
	188
	189	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
	190	IF ( wall_e_y(j,i) /= 0.0 ) THEN
	191	usvs = kappa * 0.5 * ( u(k,j,i) + u(k,j,i+1) ) &
	192	/ LOG( 0.5 * dy / z0(j,i) )
	193	usvs = usvs * ABS( usvs )
	194	dudy = wall_e_y(j,i) * usvs / km(k,j,i)
	195	ELSE
	196	dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
	197	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
	198	ENDIF
	199	dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
	200	u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
	201
	202	IF ( wall_e_x(j,i) /= 0.0 ) THEN
	203	vsus = kappa * 0.5 * ( v(k,j,i) + v(k,j+1,i) ) &
	204	/ LOG( 0.5 * dx / z0(j,i))
	205	vsus = vsus * ABS( vsus )
	206	dvdx = wall_e_x(j,i) * vsus / km(k,j,i)
	207	ELSE
	208	dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
	209	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
	210	ENDIF
	211	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
	212	dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
	213	v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
	214
	215	IF ( wall_e_x(j,i) /= 0.0 ) THEN
	216	wsus = kappa * 0.5 * ( w(k,j,i) + w(k-1,j,i) ) &
	217	/ LOG( 0.5 * dx / z0(j,i))
	218	wsus = wsus * ABS( wsus )
	219	dwdx = wall_e_x(j,i) * wsus / km(k,j,i)
	220	ELSE
	221	dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
	222	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
	223	ENDIF
	224	IF ( wall_e_y(j,i) /= 0.0 ) THEN
	225	wsvs = kappa * ( w(k,j,i) + w(k-1,j,i) ) &
	226	/ LOG( 0.5 * dy / z0(j,i))
	227	wsvs = wsvs * ABS( wsvs )
	228	dwdy = wall_e_y(j,i) * wsvs / km(k,j,i)
	229	ELSE
	230	dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
	231	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
	232	ENDIF
	233	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
	234
	235	def = 2.0 * ( dudx2 + dvdy2 + dwdz**2 ) + &
	236	dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz**2 + &
	237	dvdz*2 + 2.0 ( dvdxdudy + dwdxdudz + dwdy*dvdz )
	238
	239	IF ( def < 0.0 ) def = 0.0
	240
	241	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
	242
	243	ENDDO
	244
	245	ENDIF
	246
	247	ENDDO
	248
	249	!
	250	!-- 4 - Wird immer ausgefuehrt.
	251	!-- 'Sonderfall: Freie Atmosphaere' (wie bei 0)
	252	DO j = nys, nyn
	253
	254	IF ( ( wall_e_x(j,i) /= 0.0 ) .OR. ( wall_e_y(j,i) /= 0.0 ) ) &
	255	THEN
	256
	257	k = nzb_diff_s_outer(j,i)-1
	258
	259	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
	260	dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
	261	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
	262	dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
	263	u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
	264
	265	dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
	266	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
	267	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
	268	dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
	269	v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
	270
	271	dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
	272	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
	273	dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
	274	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
	275	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
	276
	277	def = 2.0 * ( dudx2 + dvdy2 + dwdz**2 ) + &
	278	dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz**2 + &
	279	dvdz*2 + 2.0 ( dvdxdudy + dwdxdudz + dwdy*dvdz )
	280
	281	IF ( def < 0.0 ) def = 0.0
	282
	283	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
	284
	285	ENDIF
	286
	287	ENDDO
	288
	289	!
	290	!-- Position ohne angrenzende Gebaeudewand
	291	!-- 1 - Wird immer ausgefuehrt. 'Nur Boden: u_0,v_0'
	292	DO j = nys, nyn
	293
	294	IF ( ( wall_e_x(j,i) == 0.0 ) .AND. ( wall_e_y(j,i) == 0.0 ) ) &
	295	THEN
	296
	297	k = nzb_diff_s_inner(j,i)-1
	298
	299	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
	300	dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
	301	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
	302	dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
	303	u_0(j,i) - u_0(j,i+1) ) * dd2zu(k)
	304
	305	dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
	306	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
	307	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
	308	dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
	309	v_0(j,i) - v_0(j+1,i) ) * dd2zu(k)
	310
	311	dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
	312	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
	313	dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
	314	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
	315	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
	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	ENDIF
	326
	327	ENDDO
	328
[37]	329	ELSEIF ( use_surface_fluxes ) THEN
	330
	331	DO j = nys, nyn
	332
	333	k = nzb_diff_s_outer(j,i)-1
	334
	335	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
	336	dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
	337	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
	338	dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
	339	u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
	340
	341	dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
	342	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
	343	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
	344	dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
	345	v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
	346
	347	dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
	348	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
	349	dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
	350	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
	351	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
	352
	353	def = 2.0 * ( dudx2 + dvdy2 + dwdz**2 ) + &
	354	dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz**2 + &
	355	dvdz*2 + 2.0 ( dvdxdudy + dwdxdudz + dwdy*dvdz )
	356
	357	IF ( def < 0.0 ) def = 0.0
	358
	359	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
	360
	361	ENDDO
	362
[1]	363	ENDIF
	364
	365	!
	366	!-- Calculate TKE production by buoyancy
	367	IF ( .NOT. moisture ) THEN
	368
	369	DO j = nys, nyn
	370
[19]	371	DO k = nzb_diff_s_inner(j,i), nzt_diff
[1]	372	tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / pt(k,j,i) * &
	373	( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k)
	374	ENDDO
[19]	375
[1]	376	IF ( use_surface_fluxes ) THEN
	377	k = nzb_diff_s_inner(j,i)-1
	378	tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * shf(j,i)
	379	ENDIF
	380
[19]	381	IF ( use_top_fluxes ) THEN
	382	k = nzt
	383	tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * tswst(j,i)
	384	ENDIF
	385
[1]	386	ENDDO
	387
	388	ELSE
	389
	390	DO j = nys, nyn
	391
[19]	392	DO k = nzb_diff_s_inner(j,i), nzt_diff
[1]	393
	394	IF ( .NOT. cloud_physics ) THEN
	395	k1 = 1.0 + 0.61 * q(k,j,i)
	396	k2 = 0.61 * pt(k,j,i)
	397	ELSE
	398	IF ( ql(k,j,i) == 0.0 ) THEN
	399	k1 = 1.0 + 0.61 * q(k,j,i)
	400	k2 = 0.61 * pt(k,j,i)
	401	ELSE
	402	theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
	403	temp = theta * t_d_pt(k)
	404	k1 = ( 1.0 - q(k,j,i) + 1.61 * &
	405	( q(k,j,i) - ql(k,j,i) ) * &
	406	( 1.0 + 0.622 * l_d_r / temp ) ) / &
	407	( 1.0 + 0.622 * l_d_r * l_d_cp * &
	408	( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
	409	k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
	410	ENDIF
	411	ENDIF
	412
	413	tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) * &
	414	( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + &
	415	k2 * ( q(k+1,j,i) - q(k-1,j,i) ) &
	416	) * dd2zu(k)
	417	ENDDO
	418
	419	ENDDO
	420
	421	IF ( use_surface_fluxes ) THEN
	422
	423	DO j = nys, nyn
	424
[19]	425	k = nzb_diff_s_inner(j,i)
[1]	426
	427	IF ( .NOT. cloud_physics ) THEN
	428	k1 = 1.0 + 0.61 * q(k,j,i)
	429	k2 = 0.61 * pt(k,j,i)
	430	ELSE
	431	IF ( ql(k,j,i) == 0.0 ) THEN
	432	k1 = 1.0 + 0.61 * q(k,j,i)
	433	k2 = 0.61 * pt(k,j,i)
	434	ELSE
	435	theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
	436	temp = theta * t_d_pt(k)
	437	k1 = ( 1.0 - q(k,j,i) + 1.61 * &
	438	( q(k,j,i) - ql(k,j,i) ) * &
	439	( 1.0 + 0.622 * l_d_r / temp ) ) / &
	440	( 1.0 + 0.622 * l_d_r * l_d_cp * &
	441	( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
	442	k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
	443	ENDIF
	444	ENDIF
	445
	446	tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * &
	447	( k1* shf(j,i) + k2 * qsws(j,i) )
	448	ENDDO
	449
	450	ENDIF
	451
[19]	452	IF ( use_top_fluxes ) THEN
	453
	454	DO j = nys, nyn
	455
	456	k = nzt
	457
	458	IF ( .NOT. cloud_physics ) THEN
	459	k1 = 1.0 + 0.61 * q(k,j,i)
	460	k2 = 0.61 * pt(k,j,i)
	461	ELSE
	462	IF ( ql(k,j,i) == 0.0 ) THEN
	463	k1 = 1.0 + 0.61 * q(k,j,i)
	464	k2 = 0.61 * pt(k,j,i)
	465	ELSE
	466	theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
	467	temp = theta * t_d_pt(k)
	468	k1 = ( 1.0 - q(k,j,i) + 1.61 * &
	469	( q(k,j,i) - ql(k,j,i) ) * &
	470	( 1.0 + 0.622 * l_d_r / temp ) ) / &
	471	( 1.0 + 0.622 * l_d_r * l_d_cp * &
	472	( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
	473	k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
	474	ENDIF
	475	ENDIF
	476
	477	tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * &
	478	( k1* tswst(j,i) + k2 * qswst(j,i) )
	479	ENDDO
	480
	481	ENDIF
	482
[1]	483	ENDIF
	484
	485	ENDDO
	486
	487	END SUBROUTINE production_e
	488
	489
	490	!------------------------------------------------------------------------------!
	491	! Call for grid point i,j
	492	!------------------------------------------------------------------------------!
	493	SUBROUTINE production_e_ij( i, j )
	494
	495	USE arrays_3d
	496	USE cloud_parameters
	497	USE control_parameters
	498	USE grid_variables
	499	USE indices
	500	USE statistics
	501
	502	IMPLICIT NONE
	503
	504	INTEGER :: i, j, k
	505
	506	REAL :: def, dudx, dudy, dudz, dvdx, dvdy, dvdz, dwdx, dwdy, dwdz, &
	507	k1, k2, theta, temp, usvs, vsus, wsus,wsvs
	508
	509	!
	510	!-- Calculate TKE production by shear
[19]	511	DO k = nzb_diff_s_outer(j,i), nzt
[1]	512
	513	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
	514	dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
	515	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
	516	dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
	517	u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
	518
	519	dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
	520	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
	521	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
	522	dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
	523	v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
	524
	525	dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
	526	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
	527	dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
	528	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
	529	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
	530
	531	def = 2.0 * ( dudx2 + dvdy2 + dwdz**2 ) &
	532	+ dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz2 + dvdz2 &
	533	+ 2.0 * ( dvdxdudy + dwdxdudz + dwdy*dvdz )
	534
	535	IF ( def < 0.0 ) def = 0.0
	536
	537	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
	538
	539	ENDDO
	540
[37]	541	IF ( prandtl_layer ) THEN
[1]	542
	543	IF ( ( wall_e_x(j,i) /= 0.0 ) .OR. ( wall_e_y(j,i) /= 0.0 ) ) THEN
	544	!
	545	!-- Position neben Gebaeudewand
	546	!-- 2 - Wird immer ausgefuehrt. 'Boden und Wand:
	547	!-- u_0,v_0 und Wall functions'
	548	k = nzb_diff_s_inner(j,i)-1
	549
	550	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
	551	IF ( wall_e_y(j,i) /= 0.0 ) THEN
	552	usvs = kappa * 0.5 * ( u(k,j,i) + u(k,j,i+1) ) &
	553	/ LOG( 0.5 * dy / z0(j,i) )
	554	usvs = usvs * ABS( usvs )
	555	dudy = wall_e_y(j,i) * usvs / km(k,j,i)
	556	ELSE
	557	dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
	558	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
	559	ENDIF
	560	dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
	561	u_0(j,i) - u_0(j,i+1) ) * dd2zu(k)
	562
	563	IF ( wall_e_x(j,i) /= 0.0 ) THEN
	564	vsus = kappa * 0.5 * ( v(k,j,i) + v(k,j+1,i) ) &
	565	/ LOG( 0.5 * dx / z0(j,i))
	566	vsus = vsus * ABS( vsus )
	567	dvdx = wall_e_x(j,i) * vsus / km(k,j,i)
	568	ELSE
	569	dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
	570	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
	571	ENDIF
	572	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
	573	dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
	574	v_0(j,i) - v_0(j+1,i) ) * dd2zu(k)
	575
	576	IF ( wall_e_x(j,i) /= 0.0 ) THEN
	577	wsus = kappa * 0.5 * ( w(k,j,i) + w(k-1,j,i) ) &
	578	/ LOG( 0.5 * dx / z0(j,i))
	579	wsus = wsus * ABS( wsus )
	580	dwdx = wall_e_x(j,i) * wsus / km(k,j,i)
	581	ELSE
	582	dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
	583	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
	584	ENDIF
	585	IF ( wall_e_y(j,i) /= 0.0 ) THEN
	586	wsvs = kappa * ( w(k,j,i) + w(k-1,j,i) ) &
	587	/ LOG( 0.5 * dy / z0(j,i))
	588	wsvs = wsvs * ABS( wsvs )
	589	dwdy = wall_e_y(j,i) * wsvs / km(k,j,i)
	590	ELSE
	591	dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
	592	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
	593	ENDIF
	594	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
	595
	596	def = 2.0 * ( dudx2 + dvdy2 + dwdz**2 ) + &
	597	dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz**2 + &
	598	dvdz*2 + 2.0 ( dvdxdudy + dwdxdudz + dwdy*dvdz )
	599
	600	IF ( def < 0.0 ) def = 0.0
	601
	602	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
	603
	604	!
	605	!-- 3 - Wird nur ausgefuehrt, wenn mindestens ein Niveau
	606	!-- zwischen 2 und 4 liegt, d.h. ab einer Gebaeudemindest-
	607	!-- hoehe von 2 dz. 'Nur Wand: Wall functions'
	608	DO k = nzb_diff_s_inner(j,i), nzb_diff_s_outer(j,i)-2
	609
	610	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
	611	IF ( wall_e_y(j,i) /= 0.0 ) THEN
	612	usvs = kappa * 0.5 * ( u(k,j,i) + u(k,j,i+1) ) &
	613	/ LOG( 0.5 * dy / z0(j,i) )
	614	usvs = usvs * ABS( usvs )
	615	dudy = wall_e_y(j,i) * usvs / km(k,j,i)
	616	ELSE
	617	dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
	618	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
	619	ENDIF
	620	dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
	621	u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
	622
	623	IF ( wall_e_x(j,i) /= 0.0 ) THEN
	624	vsus = kappa * 0.5 * ( v(k,j,i) + v(k,j+1,i) ) &
	625	/ LOG( 0.5 * dx / z0(j,i))
	626	vsus = vsus * ABS( vsus )
	627	dvdx = wall_e_x(j,i) * vsus / km(k,j,i)
	628	ELSE
	629	dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
	630	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
	631	ENDIF
	632	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
	633	dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
	634	v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
	635
	636	IF ( wall_e_x(j,i) /= 0.0 ) THEN
	637	wsus = kappa * 0.5 * ( w(k,j,i) + w(k-1,j,i) ) &
	638	/ LOG( 0.5 * dx / z0(j,i))
	639	wsus = wsus * ABS( wsus )
	640	dwdx = wall_e_x(j,i) * wsus / km(k,j,i)
	641	ELSE
	642	dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
	643	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
	644	ENDIF
	645	IF ( wall_e_y(j,i) /= 0.0 ) THEN
	646	wsvs = kappa * ( w(k,j,i) + w(k-1,j,i) ) &
	647	/ LOG( 0.5 * dy / z0(j,i))
	648	wsvs = wsvs * ABS( wsvs )
	649	dwdy = wall_e_y(j,i) * wsvs / km(k,j,i)
	650	ELSE
	651	dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
	652	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
	653	ENDIF
	654	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
	655
	656	def = 2.0 * ( dudx2 + dvdy2 + dwdz**2 ) + &
	657	dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz**2 + &
	658	dvdz*2 + 2.0 ( dvdxdudy + dwdxdudz + dwdy*dvdz )
	659
	660	IF ( def < 0.0 ) def = 0.0
	661
	662	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
	663
	664	ENDDO
	665
	666	!
	667	!-- 4 - Wird immer ausgefuehrt.
	668	!-- 'Sonderfall: Freie Atmosphaere' (wie bei 0)
	669	k = nzb_diff_s_outer(j,i)-1
	670
	671	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
	672	dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
	673	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
	674	dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
	675	u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
	676
	677	dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
	678	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
	679	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
	680	dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
	681	v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
	682
	683	dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
	684	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
	685	dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
	686	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
	687	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
	688
	689	def = 2.0 * ( dudx2 + dvdy2 + dwdz**2 ) + &
	690	dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz**2 + &
	691	dvdz*2 + 2.0 ( dvdxdudy + dwdxdudz + dwdy*dvdz )
	692
	693	IF ( def < 0.0 ) def = 0.0
	694
	695	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
	696
	697	ELSE
	698
	699	!
	700	!-- Position ohne angrenzende Gebaeudewand
	701	!-- 1 - Wird immer ausgefuehrt. 'Nur Boden: u_0,v_0'
	702	k = nzb_diff_s_inner(j,i)-1
	703
	704	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
	705	dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
	706	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
	707	dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
	708	u_0(j,i) - u_0(j,i+1) ) * dd2zu(k)
	709
	710	dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
	711	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
	712	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
	713	dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
	714	v_0(j,i) - v_0(j+1,i) ) * dd2zu(k)
	715
	716	dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
	717	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
	718	dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
	719	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
	720	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
	721
	722	def = 2.0 * ( dudx2 + dvdy2 + dwdz**2 ) &
	723	+ dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz2 + dvdz2 &
	724	+ 2.0 * ( dvdxdudy + dwdxdudz + dwdy*dvdz )
	725
	726	IF ( def < 0.0 ) def = 0.0
	727
	728	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
	729
	730	ENDIF
	731
[37]	732	ELSEIF ( use_surface_fluxes ) THEN
	733
	734	k = nzb_diff_s_outer(j,i)-1
	735
	736	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
	737	dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
	738	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
	739	dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
	740	u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
	741
	742	dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
	743	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
	744	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
	745	dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
	746	v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
	747
	748	dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
	749	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
	750	dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
	751	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
	752	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
	753
	754	def = 2.0 * ( dudx2 + dvdy2 + dwdz**2 ) + &
	755	dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz**2 + &
	756	dvdz*2 + 2.0 ( dvdxdudy + dwdxdudz + dwdy*dvdz )
	757
	758	IF ( def < 0.0 ) def = 0.0
	759
	760	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
	761
[1]	762	ENDIF
	763
	764	!
	765	!-- Calculate TKE production by buoyancy
	766	IF ( .NOT. moisture ) THEN
	767
[19]	768	DO k = nzb_diff_s_inner(j,i), nzt_diff
[1]	769	tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / pt(k,j,i) * &
	770	( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k)
	771	ENDDO
[19]	772
[1]	773	IF ( use_surface_fluxes ) THEN
	774	k = nzb_diff_s_inner(j,i)-1
	775	tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * shf(j,i)
	776	ENDIF
	777
[19]	778	IF ( use_top_fluxes ) THEN
	779	k = nzt
	780	tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * tswst(j,i)
	781	ENDIF
	782
[1]	783	ELSE
	784
[19]	785	DO k = nzb_diff_s_inner(j,i), nzt_diff
[1]	786
	787	IF ( .NOT. cloud_physics ) THEN
	788	k1 = 1.0 + 0.61 * q(k,j,i)
	789	k2 = 0.61 * pt(k,j,i)
	790	ELSE
	791	IF ( ql(k,j,i) == 0.0 ) THEN
	792	k1 = 1.0 + 0.61 * q(k,j,i)
	793	k2 = 0.61 * pt(k,j,i)
	794	ELSE
	795	theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
	796	temp = theta * t_d_pt(k)
	797	k1 = ( 1.0 - q(k,j,i) + 1.61 * &
	798	( q(k,j,i) - ql(k,j,i) ) * &
	799	( 1.0 + 0.622 * l_d_r / temp ) ) / &
	800	( 1.0 + 0.622 * l_d_r * l_d_cp * &
	801	( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
	802	k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
	803	ENDIF
	804	ENDIF
	805
	806	tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) * &
	807	( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + &
	808	k2 * ( q(k+1,j,i) - q(k-1,j,i) ) &
	809	) * dd2zu(k)
	810	ENDDO
[19]	811
[1]	812	IF ( use_surface_fluxes ) THEN
	813	k = nzb_diff_s_inner(j,i)-1
	814
	815	IF ( .NOT. cloud_physics ) THEN
	816	k1 = 1.0 + 0.61 * q(k,j,i)
	817	k2 = 0.61 * pt(k,j,i)
	818	ELSE
	819	IF ( ql(k,j,i) == 0.0 ) THEN
	820	k1 = 1.0 + 0.61 * q(k,j,i)
	821	k2 = 0.61 * pt(k,j,i)
	822	ELSE
	823	theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
	824	temp = theta * t_d_pt(k)
	825	k1 = ( 1.0 - q(k,j,i) + 1.61 * &
	826	( q(k,j,i) - ql(k,j,i) ) * &
	827	( 1.0 + 0.622 * l_d_r / temp ) ) / &
	828	( 1.0 + 0.622 * l_d_r * l_d_cp * &
	829	( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
	830	k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
	831	ENDIF
	832	ENDIF
	833
	834	tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * &
	835	( k1* shf(j,i) + k2 * qsws(j,i) )
	836	ENDIF
	837
[19]	838	IF ( use_top_fluxes ) THEN
	839	k = nzt
	840
	841	IF ( .NOT. cloud_physics ) THEN
	842	k1 = 1.0 + 0.61 * q(k,j,i)
	843	k2 = 0.61 * pt(k,j,i)
	844	ELSE
	845	IF ( ql(k,j,i) == 0.0 ) THEN
	846	k1 = 1.0 + 0.61 * q(k,j,i)
	847	k2 = 0.61 * pt(k,j,i)
	848	ELSE
	849	theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
	850	temp = theta * t_d_pt(k)
	851	k1 = ( 1.0 - q(k,j,i) + 1.61 * &
	852	( q(k,j,i) - ql(k,j,i) ) * &
	853	( 1.0 + 0.622 * l_d_r / temp ) ) / &
	854	( 1.0 + 0.622 * l_d_r * l_d_cp * &
	855	( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
	856	k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
	857	ENDIF
	858	ENDIF
	859
	860	tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * &
	861	( k1* tswst(j,i) + k2 * qswst(j,i) )
	862	ENDIF
	863
[1]	864	ENDIF
	865
	866	END SUBROUTINE production_e_ij
	867
	868
	869	SUBROUTINE production_e_init
	870
	871	USE arrays_3d
	872	USE control_parameters
	873	USE grid_variables
	874	USE indices
	875
	876	IMPLICIT NONE
	877
	878	INTEGER :: i, j, ku, kv
	879
[37]	880	IF ( prandtl_layer ) THEN
[1]	881
	882	IF ( first_call ) THEN
	883	ALLOCATE( u_0(nys-1:nyn+1,nxl-1:nxr+1), &
	884	v_0(nys-1:nyn+1,nxl-1:nxr+1) )
	885	first_call = .FALSE.
	886	ENDIF
	887
	888	!
	889	!-- Calculate a virtual velocity at the surface in a way that the
	890	!-- vertical velocity gradient at k = 1 (u(k+1)-u_0) matches the
	891	!-- Prandtl law (-w'u'/km). This gradient is used in the TKE shear
	892	!-- production term at k=1 (see production_e_ij).
	893	!-- The velocity gradient has to be limited in case of too small km
	894	!-- (otherwise the timestep may be significantly reduced by large
	895	!-- surface winds).
	896	!-- WARNING: the exact analytical solution would require the determination
	897	!-- of the eddy diffusivity by km = u* * kappa * zp / phi_m.
	898	!$OMP PARALLEL DO PRIVATE( ku, kv )
	899	DO i = nxl, nxr
	900	DO j = nys, nyn
	901
	902	ku = nzb_u_inner(j,i)+1
	903	kv = nzb_v_inner(j,i)+1
	904
	905	u_0(j,i) = u(ku+1,j,i) + usws(j,i) * ( zu(ku+1) - zu(ku-1) ) / &
	906	( 0.5 * ( km(ku,j,i) + km(ku,j,i-1) ) + &
	907	1.0E-20 )
	908	! ( us(j,i) * kappa * zu(1) )
	909	v_0(j,i) = v(kv+1,j,i) + vsws(j,i) * ( zu(kv+1) - zu(kv-1) ) / &
	910	( 0.5 * ( km(kv,j,i) + km(kv,j-1,i) ) + &
	911	1.0E-20 )
	912	! ( us(j,i) * kappa * zu(1) )
	913
	914	IF ( ABS( u(ku+1,j,i) - u_0(j,i) ) > &
	915	ABS( u(ku+1,j,i) - u(ku-1,j,i) ) ) u_0(j,i) = u(ku-1,j,i)
	916	IF ( ABS( v(kv+1,j,i) - v_0(j,i) ) > &
	917	ABS( v(kv+1,j,i) - v(kv-1,j,i) ) ) v_0(j,i) = v(kv-1,j,i)
	918
	919	ENDDO
	920	ENDDO
	921
	922	CALL exchange_horiz_2d( u_0 )
	923	CALL exchange_horiz_2d( v_0 )
	924
	925	ENDIF
	926
	927	END SUBROUTINE production_e_init
	928
	929	END MODULE production_e_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |