Home

Context Navigation

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

Last change on this file since 336 was 226, checked in by raasch, 16 years ago
preparations for the next release
Property svn:keywords set to `Id`
File size: 45.6 KB

Rev	Line
[1]	1	MODULE production_e_mod
	2
	3	!------------------------------------------------------------------------------!
	4	! Actual revisions:
	5	! -----------------
[226]	6	!
[110]	7	!
	8	! Former revisions:
	9	! -----------------
	10	! $Id: production_e.f90 226 2009-02-02 07:39:34Z raasch $
	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
	187	dudy = - wall_e_y(j,i) * usvs(k) / km_neutral
	188	dwdy = - wall_e_y(j,i) * wsvs(k) / km_neutral
[1]	189	ELSE
	190	dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
	191	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
[53]	192	dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
	193	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
[1]	194	ENDIF
	195
	196	IF ( wall_e_x(j,i) /= 0.0 ) THEN
[208]	197	!
	198	!-- Inconsistency removed: as the thermal stratification is
	199	!-- not taken into account for the evaluation of the wall
	200	!-- fluxes at vertical walls, the eddy viscosity km must not
	201	!-- be used for the evaluation of the velocity gradients dvdx
	202	!-- and dwdx
	203	!-- Note: The validity of the new method has not yet been
	204	!-- shown, as so far no suitable data for a validation
	205	!-- has been available
[53]	206	CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
	207	vsus, 0.0, 1.0, 0.0, 0.0 )
	208	CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
	209	wsus, 0.0, 0.0, 0.0, 1.0 )
[208]	210	km_neutral = kappa * ( vsus(k)2 + wsus(k)2 )*0.25 &
	211	0.5 * dx
	212	dvdx = - wall_e_x(j,i) * vsus(k) / km_neutral
	213	dwdx = - wall_e_x(j,i) * wsus(k) / km_neutral
[1]	214	ELSE
	215	dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
	216	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
	217	dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
	218	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
	219	ENDIF
	220
	221	def = 2.0 * ( dudx2 + dvdy2 + dwdz**2 ) + &
	222	dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz**2 + &
	223	dvdz*2 + 2.0 ( dvdxdudy + dwdxdudz + dwdy*dvdz )
	224
	225	IF ( def < 0.0 ) def = 0.0
	226
	227	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
	228
	229
	230	!
[55]	231	!-- (3) - will be executed only, if there is at least one level
	232	!-- between (2) and (4), i.e. the topography must have a
	233	!-- minimum height of 2 dz. Wall fluxes for this case have
	234	!-- already been calculated for (2).
	235	!-- 'wall only: use wall functions'
[1]	236
	237	DO k = nzb_diff_s_inner(j,i), nzb_diff_s_outer(j,i)-2
	238
	239	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
[53]	240	dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
	241	u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
	242	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
	243	dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
	244	v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
	245	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
	246
[1]	247	IF ( wall_e_y(j,i) /= 0.0 ) THEN
[208]	248	!
	249	!-- Inconsistency removed: as the thermal stratification
	250	!-- is not taken into account for the evaluation of the
	251	!-- wall fluxes at vertical walls, the eddy viscosity km
	252	!-- must not be used for the evaluation of the velocity
	253	!-- gradients dudy and dwdy
	254	!-- Note: The validity of the new method has not yet
	255	!-- been shown, as so far no suitable data for a
	256	!-- validation has been available
	257	km_neutral = kappa * ( usvs(k)**2 + &
	258	wsvs(k)2 )0.25 * 0.5 * dy
	259	dudy = - wall_e_y(j,i) * usvs(k) / km_neutral
	260	dwdy = - wall_e_y(j,i) * wsvs(k) / km_neutral
[1]	261	ELSE
	262	dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
	263	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
[53]	264	dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
	265	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
[1]	266	ENDIF
	267
	268	IF ( wall_e_x(j,i) /= 0.0 ) THEN
[208]	269	!
	270	!-- Inconsistency removed: as the thermal stratification
	271	!-- is not taken into account for the evaluation of the
	272	!-- wall fluxes at vertical walls, the eddy viscosity km
	273	!-- must not be used for the evaluation of the velocity
	274	!-- gradients dvdx and dwdx
	275	!-- Note: The validity of the new method has not yet
	276	!-- been shown, as so far no suitable data for a
	277	!-- validation has been available
	278	km_neutral = kappa * ( vsus(k)**2 + &
	279	wsus(k)2 )0.25 * 0.5 * dx
	280	dvdx = - wall_e_x(j,i) * vsus(k) / km_neutral
	281	dwdx = - wall_e_x(j,i) * wsus(k) / km_neutral
[1]	282	ELSE
	283	dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
	284	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
	285	dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
	286	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
	287	ENDIF
	288
	289	def = 2.0 * ( dudx2 + dvdy2 + dwdz**2 ) + &
	290	dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz**2 + &
	291	dvdz*2 + 2.0 ( dvdxdudy + dwdxdudz + dwdy*dvdz )
	292
	293	IF ( def < 0.0 ) def = 0.0
	294
	295	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
	296
	297	ENDDO
	298
	299	ENDIF
	300
	301	ENDDO
	302
	303	!
[55]	304	!-- (4) - will allways be executed.
	305	!-- 'special case: free atmosphere' (as for case (0))
[1]	306	DO j = nys, nyn
	307
	308	IF ( ( wall_e_x(j,i) /= 0.0 ) .OR. ( wall_e_y(j,i) /= 0.0 ) ) &
	309	THEN
	310
	311	k = nzb_diff_s_outer(j,i)-1
	312
	313	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
	314	dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
	315	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
	316	dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
	317	u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
	318
	319	dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
	320	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
	321	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
	322	dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
	323	v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
	324
	325	dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
	326	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
	327	dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
	328	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
	329	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
	330
	331	def = 2.0 * ( dudx2 + dvdy2 + dwdz**2 ) + &
	332	dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz**2 + &
	333	dvdz*2 + 2.0 ( dvdxdudy + dwdxdudz + dwdy*dvdz )
	334
	335	IF ( def < 0.0 ) def = 0.0
	336
	337	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
	338
	339	ENDIF
	340
	341	ENDDO
	342
	343	!
[55]	344	!-- Position without adjacent wall
	345	!-- (1) - will allways be executed.
	346	!-- 'bottom only: use u_0,v_0'
[1]	347	DO j = nys, nyn
	348
	349	IF ( ( wall_e_x(j,i) == 0.0 ) .AND. ( wall_e_y(j,i) == 0.0 ) ) &
	350	THEN
	351
	352	k = nzb_diff_s_inner(j,i)-1
	353
	354	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
	355	dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
	356	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
	357	dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
	358	u_0(j,i) - u_0(j,i+1) ) * dd2zu(k)
	359
	360	dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
	361	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
	362	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
	363	dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
	364	v_0(j,i) - v_0(j+1,i) ) * dd2zu(k)
	365
	366	dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
	367	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
	368	dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
	369	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
	370	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
	371
	372	def = 2.0 * ( dudx2 + dvdy2 + dwdz**2 ) + &
	373	dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz**2 + &
	374	dvdz*2 + 2.0 ( dvdxdudy + dwdxdudz + dwdy*dvdz )
	375
	376	IF ( def < 0.0 ) def = 0.0
	377
	378	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
	379
	380	ENDIF
	381
	382	ENDDO
	383
[37]	384	ELSEIF ( use_surface_fluxes ) THEN
	385
	386	DO j = nys, nyn
	387
	388	k = nzb_diff_s_outer(j,i)-1
	389
	390	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
	391	dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
	392	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
	393	dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
	394	u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
	395
	396	dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
	397	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
	398	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
	399	dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
	400	v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
	401
	402	dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
	403	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
	404	dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
	405	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
	406	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
	407
	408	def = 2.0 * ( dudx2 + dvdy2 + dwdz**2 ) + &
	409	dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz**2 + &
	410	dvdz*2 + 2.0 ( dvdxdudy + dwdxdudz + dwdy*dvdz )
	411
	412	IF ( def < 0.0 ) def = 0.0
	413
	414	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
	415
	416	ENDDO
	417
[1]	418	ENDIF
	419
	420	!
	421	!-- Calculate TKE production by buoyancy
[75]	422	IF ( .NOT. humidity ) THEN
[1]	423
[97]	424	IF ( use_reference ) THEN
[1]	425
[97]	426	IF ( ocean ) THEN
	427	!
	428	!-- So far in the ocean no special treatment of density flux in
	429	!-- the bottom and top surface layer
	430	DO j = nys, nyn
[124]	431	DO k = nzb_s_inner(j,i)+1, nzt
[97]	432	tend(k,j,i) = tend(k,j,i) + &
	433	kh(k,j,i) * g / prho_reference * &
	434	( rho(k+1,j,i)-rho(k-1,j,i) ) * dd2zu(k)
	435	ENDDO
	436	ENDDO
	437
	438	ELSE
	439
	440	DO j = nys, nyn
	441	DO k = nzb_diff_s_inner(j,i), nzt_diff
[106]	442	tend(k,j,i) = tend(k,j,i) - &
[97]	443	kh(k,j,i) * g / pt_reference * &
[57]	444	( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k)
[97]	445	ENDDO
	446
	447	IF ( use_surface_fluxes ) THEN
	448	k = nzb_diff_s_inner(j,i)-1
	449	tend(k,j,i) = tend(k,j,i) + g / pt_reference * shf(j,i)
	450	ENDIF
	451
	452	IF ( use_top_fluxes ) THEN
	453	k = nzt
	454	tend(k,j,i) = tend(k,j,i) + g / pt_reference * &
	455	tswst(j,i)
	456	ENDIF
[57]	457	ENDDO
	458
[97]	459	ENDIF
[57]	460
	461	ELSE
[1]	462
[97]	463	IF ( ocean ) THEN
	464	!
	465	!-- So far in the ocean no special treatment of density flux in
	466	!-- the bottom and top surface layer
	467	DO j = nys, nyn
[124]	468	DO k = nzb_s_inner(j,i)+1, nzt
[97]	469	tend(k,j,i) = tend(k,j,i) - &
	470	kh(k,j,i) * g / rho(k,j,i) * &
	471	( rho(k+1,j,i)-rho(k-1,j,i) ) * dd2zu(k)
	472	ENDDO
	473	ENDDO
	474
	475	ELSE
	476
	477	DO j = nys, nyn
	478	DO k = nzb_diff_s_inner(j,i), nzt_diff
	479	tend(k,j,i) = tend(k,j,i) - &
	480	kh(k,j,i) * g / pt(k,j,i) * &
[57]	481	( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k)
[97]	482	ENDDO
	483
	484	IF ( use_surface_fluxes ) THEN
	485	k = nzb_diff_s_inner(j,i)-1
	486	tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * shf(j,i)
	487	ENDIF
	488
	489	IF ( use_top_fluxes ) THEN
	490	k = nzt
	491	tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * tswst(j,i)
	492	ENDIF
[57]	493	ENDDO
[19]	494
[97]	495	ENDIF
[1]	496
[57]	497	ENDIF
	498
[1]	499	ELSE
	500
	501	DO j = nys, nyn
	502
[19]	503	DO k = nzb_diff_s_inner(j,i), nzt_diff
[1]	504
	505	IF ( .NOT. cloud_physics ) THEN
	506	k1 = 1.0 + 0.61 * q(k,j,i)
	507	k2 = 0.61 * pt(k,j,i)
	508	ELSE
	509	IF ( ql(k,j,i) == 0.0 ) THEN
	510	k1 = 1.0 + 0.61 * q(k,j,i)
	511	k2 = 0.61 * pt(k,j,i)
	512	ELSE
	513	theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
	514	temp = theta * t_d_pt(k)
	515	k1 = ( 1.0 - q(k,j,i) + 1.61 * &
	516	( q(k,j,i) - ql(k,j,i) ) * &
	517	( 1.0 + 0.622 * l_d_r / temp ) ) / &
	518	( 1.0 + 0.622 * l_d_r * l_d_cp * &
	519	( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
	520	k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
	521	ENDIF
	522	ENDIF
	523
	524	tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) * &
	525	( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + &
	526	k2 * ( q(k+1,j,i) - q(k-1,j,i) ) &
	527	) * dd2zu(k)
	528	ENDDO
	529
	530	ENDDO
	531
	532	IF ( use_surface_fluxes ) THEN
	533
	534	DO j = nys, nyn
	535
[53]	536	k = nzb_diff_s_inner(j,i)-1
[1]	537
	538	IF ( .NOT. cloud_physics ) THEN
	539	k1 = 1.0 + 0.61 * q(k,j,i)
	540	k2 = 0.61 * pt(k,j,i)
	541	ELSE
	542	IF ( ql(k,j,i) == 0.0 ) THEN
	543	k1 = 1.0 + 0.61 * q(k,j,i)
	544	k2 = 0.61 * pt(k,j,i)
	545	ELSE
	546	theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
	547	temp = theta * t_d_pt(k)
	548	k1 = ( 1.0 - q(k,j,i) + 1.61 * &
	549	( q(k,j,i) - ql(k,j,i) ) * &
	550	( 1.0 + 0.622 * l_d_r / temp ) ) / &
	551	( 1.0 + 0.622 * l_d_r * l_d_cp * &
	552	( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
	553	k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
	554	ENDIF
	555	ENDIF
	556
	557	tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * &
	558	( k1* shf(j,i) + k2 * qsws(j,i) )
	559	ENDDO
	560
	561	ENDIF
	562
[19]	563	IF ( use_top_fluxes ) THEN
	564
	565	DO j = nys, nyn
	566
	567	k = nzt
	568
	569	IF ( .NOT. cloud_physics ) THEN
	570	k1 = 1.0 + 0.61 * q(k,j,i)
	571	k2 = 0.61 * pt(k,j,i)
	572	ELSE
	573	IF ( ql(k,j,i) == 0.0 ) THEN
	574	k1 = 1.0 + 0.61 * q(k,j,i)
	575	k2 = 0.61 * pt(k,j,i)
	576	ELSE
	577	theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
	578	temp = theta * t_d_pt(k)
	579	k1 = ( 1.0 - q(k,j,i) + 1.61 * &
	580	( q(k,j,i) - ql(k,j,i) ) * &
	581	( 1.0 + 0.622 * l_d_r / temp ) ) / &
	582	( 1.0 + 0.622 * l_d_r * l_d_cp * &
	583	( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
	584	k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
	585	ENDIF
	586	ENDIF
	587
	588	tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * &
	589	( k1* tswst(j,i) + k2 * qswst(j,i) )
	590	ENDDO
	591
	592	ENDIF
	593
[1]	594	ENDIF
	595
	596	ENDDO
	597
	598	END SUBROUTINE production_e
	599
	600
	601	!------------------------------------------------------------------------------!
	602	! Call for grid point i,j
	603	!------------------------------------------------------------------------------!
	604	SUBROUTINE production_e_ij( i, j )
	605
	606	USE arrays_3d
	607	USE cloud_parameters
	608	USE control_parameters
	609	USE grid_variables
	610	USE indices
	611	USE statistics
	612
	613	IMPLICIT NONE
	614
	615	INTEGER :: i, j, k
	616
	617	REAL :: def, dudx, dudy, dudz, dvdx, dvdy, dvdz, dwdx, dwdy, dwdz, &
[208]	618	k1, k2, km_neutral, theta, temp
[1]	619
[56]	620	REAL, DIMENSION(nzb:nzt+1) :: usvs, vsus, wsus, wsvs
[53]	621
[1]	622	!
	623	!-- Calculate TKE production by shear
[19]	624	DO k = nzb_diff_s_outer(j,i), nzt
[1]	625
	626	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
	627	dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
	628	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
	629	dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
	630	u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
	631
	632	dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
	633	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
	634	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
	635	dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
	636	v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
	637
	638	dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
	639	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
	640	dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
	641	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
	642	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
	643
	644	def = 2.0 * ( dudx2 + dvdy2 + dwdz**2 ) &
	645	+ dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz2 + dvdz2 &
	646	+ 2.0 * ( dvdxdudy + dwdxdudz + dwdy*dvdz )
	647
	648	IF ( def < 0.0 ) def = 0.0
	649
	650	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
	651
	652	ENDDO
	653
[37]	654	IF ( prandtl_layer ) THEN
[1]	655
	656	IF ( ( wall_e_x(j,i) /= 0.0 ) .OR. ( wall_e_y(j,i) /= 0.0 ) ) THEN
[55]	657
[1]	658	!
[55]	659	!-- Position beneath wall
	660	!-- (2) - Will allways be executed.
	661	!-- 'bottom and wall: use u_0,v_0 and wall functions'
[1]	662	k = nzb_diff_s_inner(j,i)-1
	663
	664	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
[53]	665	dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
	666	u_0(j,i) - u_0(j,i+1) ) * dd2zu(k)
	667	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
	668	dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
	669	v_0(j,i) - v_0(j+1,i) ) * dd2zu(k)
	670	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
	671
[1]	672	IF ( wall_e_y(j,i) /= 0.0 ) THEN
[208]	673	!
	674	!-- Inconsistency removed: as the thermal stratification
	675	!-- is not taken into account for the evaluation of the
	676	!-- wall fluxes at vertical walls, the eddy viscosity km
	677	!-- must not be used for the evaluation of the velocity
	678	!-- gradients dudy and dwdy
	679	!-- Note: The validity of the new method has not yet
	680	!-- been shown, as so far no suitable data for a
	681	!-- validation has been available
[53]	682	CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
	683	usvs, 1.0, 0.0, 0.0, 0.0 )
	684	CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
	685	wsvs, 0.0, 0.0, 1.0, 0.0 )
[208]	686	km_neutral = kappa * ( usvs(k)2 + wsvs(k)2 )*0.25 &
	687	0.5 * dy
	688	dudy = - wall_e_y(j,i) * usvs(k) / km_neutral
	689	dwdy = - wall_e_y(j,i) * wsvs(k) / km_neutral
[1]	690	ELSE
	691	dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
	692	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
[53]	693	dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
	694	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
[1]	695	ENDIF
	696
	697	IF ( wall_e_x(j,i) /= 0.0 ) THEN
[208]	698	!
	699	!-- Inconsistency removed: as the thermal stratification
	700	!-- is not taken into account for the evaluation of the
	701	!-- wall fluxes at vertical walls, the eddy viscosity km
	702	!-- must not be used for the evaluation of the velocity
	703	!-- gradients dvdx and dwdx
	704	!-- Note: The validity of the new method has not yet
	705	!-- been shown, as so far no suitable data for a
	706	!-- validation has been available
[53]	707	CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
	708	vsus, 0.0, 1.0, 0.0, 0.0 )
	709	CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
	710	wsus, 0.0, 0.0, 0.0, 1.0 )
[208]	711	km_neutral = kappa * ( vsus(k)2 + wsus(k)2 )*0.25 &
	712	0.5 * dx
	713	dvdx = - wall_e_x(j,i) * vsus(k) / km_neutral
	714	dwdx = - wall_e_x(j,i) * wsus(k) / km_neutral
[1]	715	ELSE
	716	dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
	717	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
	718	dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
	719	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
	720	ENDIF
	721
	722	def = 2.0 * ( dudx2 + dvdy2 + dwdz**2 ) + &
	723	dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz**2 + &
	724	dvdz*2 + 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	!
[55]	731	!-- (3) - will be executed only, if there is at least one level
	732	!-- between (2) and (4), i.e. the topography must have a
	733	!-- minimum height of 2 dz. Wall fluxes for this case have
	734	!-- already been calculated for (2).
	735	!-- 'wall only: use wall functions'
[1]	736	DO k = nzb_diff_s_inner(j,i), nzb_diff_s_outer(j,i)-2
	737
	738	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
[53]	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	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
	742	dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
	743	v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
	744	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
	745
[1]	746	IF ( wall_e_y(j,i) /= 0.0 ) THEN
[208]	747	!
	748	!-- Inconsistency removed: as the thermal stratification
	749	!-- is not taken into account for the evaluation of the
	750	!-- wall fluxes at vertical walls, the eddy viscosity km
	751	!-- must not be used for the evaluation of the velocity
	752	!-- gradients dudy and dwdy
	753	!-- Note: The validity of the new method has not yet
	754	!-- been shown, as so far no suitable data for a
	755	!-- validation has been available
	756	km_neutral = kappa * ( usvs(k)**2 + &
	757	wsvs(k)2 )0.25 * 0.5 * dy
	758	dudy = - wall_e_y(j,i) * usvs(k) / km_neutral
	759	dwdy = - wall_e_y(j,i) * wsvs(k) / km_neutral
[1]	760	ELSE
	761	dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
	762	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
[53]	763	dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
	764	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
[1]	765	ENDIF
	766
	767	IF ( wall_e_x(j,i) /= 0.0 ) THEN
[208]	768	!
	769	!-- Inconsistency removed: as the thermal stratification
	770	!-- is not taken into account for the evaluation of the
	771	!-- wall fluxes at vertical walls, the eddy viscosity km
	772	!-- must not be used for the evaluation of the velocity
	773	!-- gradients dvdx and dwdx
	774	!-- Note: The validity of the new method has not yet
	775	!-- been shown, as so far no suitable data for a
	776	!-- validation has been available
	777	km_neutral = kappa * ( vsus(k)**2 + &
	778	wsus(k)2 )0.25 * 0.5 * dx
	779	dvdx = - wall_e_x(j,i) * vsus(k) / km_neutral
	780	dwdx = - wall_e_x(j,i) * wsus(k) / km_neutral
[1]	781	ELSE
	782	dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
	783	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
	784	dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
	785	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
	786	ENDIF
	787
	788	def = 2.0 * ( dudx2 + dvdy2 + dwdz**2 ) + &
	789	dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz**2 + &
	790	dvdz*2 + 2.0 ( dvdxdudy + dwdxdudz + dwdy*dvdz )
	791
	792	IF ( def < 0.0 ) def = 0.0
	793
	794	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
	795
	796	ENDDO
	797
	798	!
[55]	799	!-- (4) - will allways be executed.
	800	!-- 'special case: free atmosphere' (as for case (0))
[1]	801	k = nzb_diff_s_outer(j,i)-1
	802
	803	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
	804	dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
	805	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
	806	dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
	807	u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
	808
	809	dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
	810	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
	811	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
	812	dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
	813	v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
	814
	815	dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
	816	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
	817	dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
	818	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
	819	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
	820
	821	def = 2.0 * ( dudx2 + dvdy2 + dwdz**2 ) + &
	822	dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz**2 + &
	823	dvdz*2 + 2.0 ( dvdxdudy + dwdxdudz + dwdy*dvdz )
	824
	825	IF ( def < 0.0 ) def = 0.0
	826
	827	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
	828
	829	ELSE
	830
	831	!
[55]	832	!-- Position without adjacent wall
	833	!-- (1) - will allways be executed.
	834	!-- 'bottom only: use u_0,v_0'
[1]	835	k = nzb_diff_s_inner(j,i)-1
	836
	837	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
	838	dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
	839	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
	840	dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
	841	u_0(j,i) - u_0(j,i+1) ) * dd2zu(k)
	842
	843	dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
	844	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
	845	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
	846	dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
	847	v_0(j,i) - v_0(j+1,i) ) * dd2zu(k)
	848
	849	dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
	850	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
	851	dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
	852	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
	853	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
	854
	855	def = 2.0 * ( dudx2 + dvdy2 + dwdz**2 ) &
	856	+ dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz2 + dvdz2 &
	857	+ 2.0 * ( dvdxdudy + dwdxdudz + dwdy*dvdz )
	858
	859	IF ( def < 0.0 ) def = 0.0
	860
	861	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
	862
	863	ENDIF
	864
[37]	865	ELSEIF ( use_surface_fluxes ) THEN
	866
	867	k = nzb_diff_s_outer(j,i)-1
	868
	869	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
	870	dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
	871	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
	872	dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
	873	u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
	874
	875	dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
	876	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
	877	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
	878	dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
	879	v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
	880
	881	dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
	882	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
	883	dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
	884	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
	885	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
	886
	887	def = 2.0 * ( dudx2 + dvdy2 + dwdz**2 ) + &
	888	dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz**2 + &
	889	dvdz*2 + 2.0 ( dvdxdudy + dwdxdudz + dwdy*dvdz )
	890
	891	IF ( def < 0.0 ) def = 0.0
	892
	893	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
	894
[1]	895	ENDIF
	896
	897	!
	898	!-- Calculate TKE production by buoyancy
[75]	899	IF ( .NOT. humidity ) THEN
[1]	900
[97]	901	IF ( use_reference ) THEN
[19]	902
[97]	903	IF ( ocean ) THEN
	904	!
	905	!-- So far in the ocean no special treatment of density flux in the
	906	!-- bottom and top surface layer
[108]	907	DO k = nzb_s_inner(j,i)+1, nzt
[97]	908	tend(k,j,i) = tend(k,j,i) + kh(k,j,i) * g / prho_reference * &
	909	( rho(k+1,j,i) - rho(k-1,j,i) ) * dd2zu(k)
	910	ENDDO
	911
	912	ELSE
	913
	914	DO k = nzb_diff_s_inner(j,i), nzt_diff
	915	tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / pt_reference * &
[57]	916	( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k)
[97]	917	ENDDO
[1]	918
[97]	919	IF ( use_surface_fluxes ) THEN
	920	k = nzb_diff_s_inner(j,i)-1
	921	tend(k,j,i) = tend(k,j,i) + g / pt_reference * shf(j,i)
	922	ENDIF
[19]	923
[97]	924	IF ( use_top_fluxes ) THEN
	925	k = nzt
	926	tend(k,j,i) = tend(k,j,i) + g / pt_reference * tswst(j,i)
	927	ENDIF
	928
[57]	929	ENDIF
	930
	931	ELSE
	932
[97]	933	IF ( ocean ) THEN
	934	!
	935	!-- So far in the ocean no special treatment of density flux in the
	936	!-- bottom and top surface layer
[124]	937	DO k = nzb_s_inner(j,i)+1, nzt
[97]	938	tend(k,j,i) = tend(k,j,i) + kh(k,j,i) * g / rho(k,j,i) * &
	939	( rho(k+1,j,i) - rho(k-1,j,i) ) * dd2zu(k)
	940	ENDDO
	941
	942	ELSE
	943
	944	DO k = nzb_diff_s_inner(j,i), nzt_diff
	945	tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / pt(k,j,i) * &
[57]	946	( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k)
[97]	947	ENDDO
[57]	948
[97]	949	IF ( use_surface_fluxes ) THEN
	950	k = nzb_diff_s_inner(j,i)-1
	951	tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * shf(j,i)
	952	ENDIF
[57]	953
[97]	954	IF ( use_top_fluxes ) THEN
	955	k = nzt
	956	tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * tswst(j,i)
	957	ENDIF
	958
[57]	959	ENDIF
	960
[97]	961	ENDIF
[57]	962
[1]	963	ELSE
	964
[19]	965	DO k = nzb_diff_s_inner(j,i), nzt_diff
[1]	966
	967	IF ( .NOT. cloud_physics ) THEN
	968	k1 = 1.0 + 0.61 * q(k,j,i)
	969	k2 = 0.61 * pt(k,j,i)
	970	ELSE
	971	IF ( ql(k,j,i) == 0.0 ) THEN
	972	k1 = 1.0 + 0.61 * q(k,j,i)
	973	k2 = 0.61 * pt(k,j,i)
	974	ELSE
	975	theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
	976	temp = theta * t_d_pt(k)
	977	k1 = ( 1.0 - q(k,j,i) + 1.61 * &
	978	( q(k,j,i) - ql(k,j,i) ) * &
	979	( 1.0 + 0.622 * l_d_r / temp ) ) / &
	980	( 1.0 + 0.622 * l_d_r * l_d_cp * &
	981	( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
	982	k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
	983	ENDIF
	984	ENDIF
	985
	986	tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) * &
	987	( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + &
	988	k2 * ( q(k+1,j,i) - q(k-1,j,i) ) &
	989	) * dd2zu(k)
	990	ENDDO
[19]	991
[1]	992	IF ( use_surface_fluxes ) THEN
	993	k = nzb_diff_s_inner(j,i)-1
	994
	995	IF ( .NOT. cloud_physics ) THEN
	996	k1 = 1.0 + 0.61 * q(k,j,i)
	997	k2 = 0.61 * pt(k,j,i)
	998	ELSE
	999	IF ( ql(k,j,i) == 0.0 ) THEN
	1000	k1 = 1.0 + 0.61 * q(k,j,i)
	1001	k2 = 0.61 * pt(k,j,i)
	1002	ELSE
	1003	theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
	1004	temp = theta * t_d_pt(k)
	1005	k1 = ( 1.0 - q(k,j,i) + 1.61 * &
	1006	( q(k,j,i) - ql(k,j,i) ) * &
	1007	( 1.0 + 0.622 * l_d_r / temp ) ) / &
	1008	( 1.0 + 0.622 * l_d_r * l_d_cp * &
	1009	( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
	1010	k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
	1011	ENDIF
	1012	ENDIF
	1013
	1014	tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * &
	1015	( k1* shf(j,i) + k2 * qsws(j,i) )
	1016	ENDIF
	1017
[19]	1018	IF ( use_top_fluxes ) THEN
	1019	k = nzt
	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) + g / vpt(k,j,i) * &
	1041	( k1* tswst(j,i) + k2 * qswst(j,i) )
	1042	ENDIF
	1043
[1]	1044	ENDIF
	1045
	1046	END SUBROUTINE production_e_ij
	1047
	1048
	1049	SUBROUTINE production_e_init
	1050
	1051	USE arrays_3d
	1052	USE control_parameters
	1053	USE grid_variables
	1054	USE indices
	1055
	1056	IMPLICIT NONE
	1057
	1058	INTEGER :: i, j, ku, kv
	1059
[37]	1060	IF ( prandtl_layer ) THEN
[1]	1061
	1062	IF ( first_call ) THEN
	1063	ALLOCATE( u_0(nys-1:nyn+1,nxl-1:nxr+1), &
	1064	v_0(nys-1:nyn+1,nxl-1:nxr+1) )
	1065	first_call = .FALSE.
	1066	ENDIF
	1067
	1068	!
	1069	!-- Calculate a virtual velocity at the surface in a way that the
	1070	!-- vertical velocity gradient at k = 1 (u(k+1)-u_0) matches the
	1071	!-- Prandtl law (-w'u'/km). This gradient is used in the TKE shear
	1072	!-- production term at k=1 (see production_e_ij).
	1073	!-- The velocity gradient has to be limited in case of too small km
	1074	!-- (otherwise the timestep may be significantly reduced by large
	1075	!-- surface winds).
[106]	1076	!-- Upper bounds are nxr+1 and nyn+1 because otherwise these values are
	1077	!-- not available in case of non-cyclic boundary conditions.
[1]	1078	!-- WARNING: the exact analytical solution would require the determination
	1079	!-- of the eddy diffusivity by km = u* * kappa * zp / phi_m.
	1080	!$OMP PARALLEL DO PRIVATE( ku, kv )
[106]	1081	DO i = nxl, nxr+1
	1082	DO j = nys, nyn+1
[1]	1083
	1084	ku = nzb_u_inner(j,i)+1
	1085	kv = nzb_v_inner(j,i)+1
	1086
	1087	u_0(j,i) = u(ku+1,j,i) + usws(j,i) * ( zu(ku+1) - zu(ku-1) ) / &
	1088	( 0.5 * ( km(ku,j,i) + km(ku,j,i-1) ) + &
	1089	1.0E-20 )
	1090	! ( us(j,i) * kappa * zu(1) )
	1091	v_0(j,i) = v(kv+1,j,i) + vsws(j,i) * ( zu(kv+1) - zu(kv-1) ) / &
	1092	( 0.5 * ( km(kv,j,i) + km(kv,j-1,i) ) + &
	1093	1.0E-20 )
	1094	! ( us(j,i) * kappa * zu(1) )
	1095
	1096	IF ( ABS( u(ku+1,j,i) - u_0(j,i) ) > &
	1097	ABS( u(ku+1,j,i) - u(ku-1,j,i) ) ) u_0(j,i) = u(ku-1,j,i)
	1098	IF ( ABS( v(kv+1,j,i) - v_0(j,i) ) > &
	1099	ABS( v(kv+1,j,i) - v(kv-1,j,i) ) ) v_0(j,i) = v(kv-1,j,i)
	1100
	1101	ENDDO
	1102	ENDDO
	1103
	1104	CALL exchange_horiz_2d( u_0 )
	1105	CALL exchange_horiz_2d( v_0 )
	1106
	1107	ENDIF
	1108
	1109	END SUBROUTINE production_e_init
	1110
	1111	END MODULE production_e_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |