Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |