Home

Context Navigation

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

Last change on this file since 364 was 364, checked in by letzel, 15 years ago
Bugfix: avoid zero division by km_neutral (production_e)
Property svn:keywords set to `Id`
File size: 47.1 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |