Home

Context Navigation

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

Last change on this file since 108 was 108, checked in by letzel, 17 years ago

Improved coupler: evaporation - salinity-flux coupling for humidity = .T.,

avoid MPI hangs when coupled runs terminate, add DOC/app/chapter_3.8;

Optional calculation of km and kh from initial TKE e_init;
Default initialization of km,kh = 0.00001 for ocean = .T.;
Allow data_output_pr= q, wq, w"q", w*q* for humidity = .T.;
Bugfix: Rayleigh damping for ocean fixed.

Property svn:keywords set to Id

File size: 40.5 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |