Home

Context Navigation

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

Last change on this file since 847 was 760, checked in by raasch, 13 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 47.6 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |