Home

Context Navigation

production_e.f90 @ 90

Last change on this file since 90 was 77, checked in by raasch, 18 years ago

New:
---

particle reflection from vertical walls implemented, particle SGS model adjusted to walls

Wall functions for vertical walls now include diabatic conditions. New subroutines wall_fluxes, wall_fluxes_e. New 4D-array rif_wall.

new d3par-parameter netcdf_64bit_3d to switch on 64bit offset only for 3D files

new d3par-parameter dt_max to define the maximum value for the allowed timestep

new inipar-parameter loop_optimization to control the loop optimization method

new inipar-parameter pt_refrence. If given, this value is used as the reference that in buoyancy terms (otherwise, the instantaneous horizontally averaged temperature is used).

new user interface user_advec_particles

new initializing action "by_user" calls user_init_3d_model and allows the initial setting of all 3d arrays

topography height informations are stored on arrays zu_s_inner and zw_w_inner and output to the 2d/3d NetCDF files

samples added to the user interface which show how to add user-define time series quantities.

calculation/output of precipitation amount, precipitation rate and z0 (by setting "pra*", "prr*", "z0*" with data_output). The time interval on which the precipitation amount is defined is set by new d3par-parameter precipitation_amount_interval

unit 9 opened for debug output (file DEBUG_<pe#>)

Makefile, advec_particles, average_3d_data, buoyancy, calc_precipitation, check_open, check_parameters, data_output_2d, diffusion_e, diffusion_u, diffusion_v, diffusion_w, diffusivities, header, impact_of_latent_heat, init_particles, init_3d_model, modules, netcdf, parin, production_e, read_var_list, read_3d_binary, sum_up_3d_data, user_interface, write_var_list, write_3d_binary

New: wall_fluxes

Changed:

General revision of non-cyclic horizontal boundary conditions: radiation boundary conditions are now used instead of Neumann conditions at the outflow (calculation needs velocity values for t-dt, which are stored on new arrays u_m_l, u_m_r, etc.), calculation of mean outflow is not needed any more, volume flow control is added for the outflow boundary (currently only for the north boundary!!), additional gridpoints along x and y (uxrp, vynp) are not needed any more, routine "boundary_conds" now operates on timelevel t+dt and is not split in two parts (main, uvw_outflow) any more, Neumann boundary conditions at inflow/outflow in case of non-cyclic boundary conditions for all 2d-arrays that are handled by exchange_horiz_2d

The FFT-method for solving the Poisson-equation is now working with Neumann boundary conditions both at the bottom and the top. This requires adjustments of the tridiagonal coefficients and subtracting the horizontally averaged mean from the vertical velocity field.

+age_m in particle_type

Particles-package is now part of the default code ("-p particles" is not needed any more).

Move call of user_actions( 'after_integration' ) below increment of times
and counters. user_actions is now called for each statistic region and has as an argument the number of the respective region (sr)

d3par-parameter data_output_ts removed. Timeseries output for "profil" removed. Timeseries are now switched on by dt_dots. Timeseries data is collected in flow_statistics.

Initial velocities at nzb+1 are regarded for volume flow control in case they have been set zero before (to avoid small timesteps); see new internal parameters u/v_nzb_p1_for_vfc.

q is not allowed to become negative (prognostic_equations).

poisfft_init is only called if fft-solver is switched on (init_pegrid).

d3par-parameter moisture renamed to humidity.

Subversion global revision number is read from mrun and added to the run description header and to the run control (_rc) file.

vtk directives removed from main program.

The uitility routine interpret_config reads PALM environment variables from NAMELIST instead using the system call GETENV.

advec_u_pw, advec_u_up, advec_v_pw, advec_v_up, asselin_filter, check_parameters, coriolis, data_output_dvrp, data_output_ptseries, data_output_ts, data_output_2d, data_output_3d, diffusion_u, diffusion_v, exchange_horiz, exchange_horiz_2d, flow_statistics, header, init_grid, init_particles, init_pegrid, init_rankine, init_pt_anomaly, init_1d_model, init_3d_model, modules, palm, package_parin, parin, poisfft, poismg, prandtl_fluxes, pres, production_e, prognostic_equations, read_var_list, read_3d_binary, sor, swap_timelevel, time_integration, write_var_list, write_3d_binary

Errors:

Bugfix: preset of tendencies te_em, te_um, te_vm in init_1d_model

Bugfix in sample for reading user defined data from restart file (user_init)

Bugfix in setting diffusivities for cases with the outflow damping layer extending over more than one subdomain (init_3d_model)

Check for possible negative humidities in the initial humidity profile.

in Makefile, default suffixes removed from the suffix list to avoid calling of m2c in
# case of .mod files

Makefile
check_parameters, init_1d_model, init_3d_model, user_interface

Property svn:keywords set to Id

File size: 37.5 KB

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

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

Context Navigation

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

Download in other formats: