Home

Context Navigation

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

Last change on this file since 196 was 187, checked in by letzel, 16 years ago
new: descriptions of plant canopy model and prandtl layer (trunk/DOC/misc) changed: more consistent flux definitions; modification of the integrated version of the profile function for momentum for unstable stratification (wall_fluxes, production_e) bugfix: change definition of us_wall from 1D to 2D (prandtl_fluxes, wall_fluxes)
Property svn:keywords set to `Id`
File size: 40.6 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |