Home

Context Navigation

production_e.f90 @ 150

Last change on this file since 150 was 139, checked in by raasch, 17 years ago

New:
---

Plant canopy model of Watanabe (2004,BLM 112,307-341) added.
It can be switched on by the inipar parameter plant_canopy.
The inipar parameter canopy_mode can be used to prescribe a
plant canopy type. The default case is a homogeneous plant
canopy. Heterogeneous distributions of the leaf area
density and the canopy drag coefficient can be defined in the
new routine user_init_plant_canopy (user_interface).
The inipar parameters lad_surface, lad_vertical_gradient and
lad_vertical_gradient_level can be used in order to
prescribe the vertical profile of leaf area density. The
inipar parameter drag_coefficient determines the canopy
drag coefficient.
Finally, the inipar parameter pch_index determines the
index of the upper boundary of the plant canopy.

Allow new case bc_uv_t = 'dirichlet_0' for channel flow.

For unknown variables (CASE DEFAULT) call new subroutine user_data_output_dvrp

Pressure boundary conditions for vertical walls added to the multigrid solver.
They are applied using new wall flag arrays (wall_flags_..) which are defined
for each grid level. New argument gls added to routine user_init_grid
(user_interface).

Frequence of sorting particles can be controlled with new particles_par
parameter dt_sort_particles. Sorting is moved from the SGS timestep loop in
advec_particles after the end of this loop.

advec_particles, check_parameters, data_output_dvrp, header, init_3d_model, init_grid, init_particles, init_pegrid, modules, package_parin, parin, plant_canopy_model, read_var_list, read_3d_binary, user_interface, write_var_list, write_3d_binary

Changed:

Redefine initial nzb_local as the actual total size of topography (later the
extent of topography in nzb_local is reduced by 1dx at the E topography walls
and by 1dy at the N topography walls to form the basis for nzb_s_inner);
for consistency redefine 'single_building' case.

Vertical profiles now based on nzb_s_inner; they are divided by
ngp_2dh_s_inner (scalars, procucts of scalars) and ngp_2dh (staggered velocity
components and their products, procucts of scalars and velocity components),
respectively.

Allow two instead of one digit to specify isosurface and slicer variables.

Status of 3D-volume NetCDF data file only depends on switch netcdf_64bit_3d (check_open)

prognostic_equations include the respective wall_*flux in the parameter list of
calls of diffusion_s. Same as before, only the values of wall_heatflux(0:4)
can be assigned. At present, wall_humidityflux, wall_qflux, wall_salinityflux,
and wall_scalarflux are kept zero. diffusion_s uses the respective wall_*flux
instead of wall_heatflux. This update serves two purposes:

it avoids errors in calculations with humidity/scalar/salinity and prescribed

non-zero wall_heatflux,

it prepares PALM for a possible assignment of wall fluxes of

humidity/scalar/salinity in a future release.

buoyancy, check_open, data_output_dvrp, diffusion_s, diffusivities, flow_statistics, header, init_3d_model, init_dvrp, init_grid, modules, prognostic_equations

Errors:

Bugfix: summation of sums_l_l in diffusivities.

Several bugfixes in the ocean part: Initial density rho is calculated
(init_ocean). Error in initializing u_init and v_init removed
(check_parameters). Calculation of density flux now starts from
nzb+1 (production_e).

Bugfix: pleft/pright changed to pnorth/psouth in sendrecv of particle tail
numbers along y, small bugfixes in the SGS part (advec_particles)

Bugfix: model_string needed a default value (combine_plot_fields)

Bugfix: wavenumber calculation for even nx in routines maketri (poisfft)

Bugfix: assignment of fluxes at walls

Bugfix: absolute value of f must be used when calculating the Blackadar mixing length (init_1d_model)

advec_particles, check_parameters, combine_plot_fields, diffusion_s, diffusivities, init_ocean, init_1d_model, poisfft, production_e

Property svn:keywords set to Id

File size: 40.4 KB

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

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

Context Navigation

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

Download in other formats: