Home

Context Navigation

source: palm/trunk/SOURCE/diffusion_e.f90 @ 1332

Last change on this file since 1332 was 1321, checked in by raasch, 11 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 24.1 KB

Rev	Line
[1]	1	MODULE diffusion_e_mod
	2
[1036]	3	!--------------------------------------------------------------------------------!
	4	! This file is part of PALM.
	5	!
	6	! PALM is free software: you can redistribute it and/or modify it under the terms
	7	! of the GNU General Public License as published by the Free Software Foundation,
	8	! either version 3 of the License, or (at your option) any later version.
	9	!
	10	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
	11	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
	12	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
	13	!
	14	! You should have received a copy of the GNU General Public License along with
	15	! PALM. If not, see <http://www.gnu.org/licenses/>.
	16	!
[1310]	17	! Copyright 1997-2014 Leibniz Universitaet Hannover
[1036]	18	!--------------------------------------------------------------------------------!
	19	!
[484]	20	! Current revisions:
[1]	21	! -----------------
[1321]	22	!
	23	!
	24	! Former revisions:
	25	! -----------------
	26	! $Id: diffusion_e.f90 1321 2014-03-20 09:40:40Z suehring $
	27	!
	28	! 1320 2014-03-20 08:40:49Z raasch
[1320]	29	! ONLY-attribute added to USE-statements,
	30	! kind-parameters added to all INTEGER and REAL declaration statements,
	31	! kinds are defined in new module kinds,
	32	! revision history before 2012 removed,
	33	! comment fields (!:) to be used for variable explanations added to
	34	! all variable declaration statements
[98]	35	!
[1258]	36	! 1257 2013-11-08 15:18:40Z raasch
	37	! openacc loop and loop vector clauses removed
	38	!
[1182]	39	! 1179 2013-06-14 05:57:58Z raasch
	40	! use_reference renamed use_single_reference_value
	41	!
[1172]	42	! 1171 2013-05-30 11:27:45Z raasch
	43	! use_reference-case activated in accelerator version
	44	!
[1132]	45	! 1128 2013-04-12 06:19:32Z raasch
	46	! loop index bounds in accelerator version replaced by i_left, i_right, j_south,
	47	! j_north
	48	!
[1066]	49	! 1065 2012-11-22 17:42:36Z hoffmann
	50	! Enabled the claculation of diss in case of turbulence = .TRUE. (parameterized
	51	! effects of turbulence on autoconversion and accretion in two-moments cloud
	52	! physics scheme).
	53	!
[1037]	54	! 1036 2012-10-22 13:43:42Z raasch
	55	! code put under GPL (PALM 3.9)
	56	!
[1017]	57	! 1015 2012-09-27 09:23:24Z raasch
	58	! accelerator version (*_acc) added,
	59	! adjustment of mixing length to the Prandtl mixing length at first grid point
	60	! above ground removed
	61	!
[1011]	62	! 1010 2012-09-20 07:59:54Z raasch
	63	! cpp switch __nopointer added for pointer free version
	64	!
[1002]	65	! 1001 2012-09-13 14:08:46Z raasch
	66	! most arrays comunicated by module instead of parameter list
	67	!
[826]	68	! 825 2012-02-19 03:03:44Z raasch
	69	! wang_collision_kernel renamed wang_kernel
	70	!
[1]	71	! Revision 1.1 1997/09/19 07:40:24 raasch
	72	! Initial revision
	73	!
	74	!
	75	! Description:
	76	! ------------
	77	! Diffusion- and dissipation terms for the TKE
	78	!------------------------------------------------------------------------------!
	79
	80	PRIVATE
[1015]	81	PUBLIC diffusion_e, diffusion_e_acc
[1]	82
	83
	84	INTERFACE diffusion_e
	85	MODULE PROCEDURE diffusion_e
	86	MODULE PROCEDURE diffusion_e_ij
	87	END INTERFACE diffusion_e
	88
[1015]	89	INTERFACE diffusion_e_acc
	90	MODULE PROCEDURE diffusion_e_acc
	91	END INTERFACE diffusion_e_acc
	92
[1]	93	CONTAINS
	94
	95
	96	!------------------------------------------------------------------------------!
	97	! Call for all grid points
	98	!------------------------------------------------------------------------------!
[1001]	99	SUBROUTINE diffusion_e( var, var_reference )
[1]	100
[1320]	101	USE arrays_3d, &
	102	ONLY: dd2zu, ddzu, ddzw, diss, e, km, l_grid, tend, zu, zw
	103
	104	USE control_parameters, &
	105	ONLY: atmos_ocean_sign, g, turbulence, use_single_reference_value, &
	106	wall_adjustment, wall_adjustment_factor
	107
	108	USE grid_variables, &
	109	ONLY: ddx2, ddy2
	110
	111	USE indices, &
	112	ONLY: nxl, nxr, nyn, nys, nzb, nzb_s_inner, nzt
	113
	114	USE kinds
	115
	116	USE particle_attributes, &
	117	ONLY: use_sgs_for_particles, wang_kernel
[1]	118
	119	IMPLICIT NONE
	120
[1320]	121	INTEGER(iwp) :: i !:
	122	INTEGER(iwp) :: j !:
	123	INTEGER(iwp) :: k !:
	124	REAL(wp) :: dvar_dz !:
	125	REAL(wp) :: l_stable !:
	126	REAL(wp) :: var_reference !:
[1001]	127
[1010]	128	#if defined( __nopointer )
[1320]	129	REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: var !:
[1010]	130	#else
[1320]	131	REAL(wp), DIMENSION(:,:,:), POINTER :: var !:
[1010]	132	#endif
[1320]	133	REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: dissipation !:
	134	REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: l !:
	135	REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: ll !:
[1]	136
	137
	138	!
[65]	139	!-- This if clause must be outside the k-loop because otherwise
	140	!-- runtime errors occur with -C hopt on NEC
[1179]	141	IF ( use_single_reference_value ) THEN
[65]	142
	143	DO i = nxl, nxr
	144	DO j = nys, nyn
	145	DO k = nzb_s_inner(j,i)+1, nzt
[1]	146	!
[65]	147	!-- Calculate the mixing length (for dissipation)
[97]	148	dvar_dz = atmos_ocean_sign * &
	149	( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
	150	IF ( dvar_dz > 0.0 ) THEN
[57]	151	l_stable = 0.76 * SQRT( e(k,j,i) ) / &
[97]	152	SQRT( g / var_reference * dvar_dz ) + 1E-5
[57]	153	ELSE
[65]	154	l_stable = l_grid(k)
[57]	155	ENDIF
[1]	156	!
[65]	157	!-- Adjustment of the mixing length
	158	IF ( wall_adjustment ) THEN
[94]	159	l(k,j) = MIN( wall_adjustment_factor * &
	160	( zu(k) - zw(nzb_s_inner(j,i)) ), &
	161	l_grid(k), l_stable )
	162	ll(k,j) = MIN( wall_adjustment_factor * &
	163	( zu(k) - zw(nzb_s_inner(j,i)) ), &
	164	l_grid(k) )
[65]	165	ELSE
	166	l(k,j) = MIN( l_grid(k), l_stable )
	167	ll(k,j) = l_grid(k)
	168	ENDIF
[1]	169
[65]	170	ENDDO
[1]	171	ENDDO
[65]	172
[1]	173	!
[65]	174	!-- Calculate the tendency terms
	175	DO j = nys, nyn
	176	DO k = nzb_s_inner(j,i)+1, nzt
[1]	177
[65]	178	dissipation(k,j) = ( 0.19 + 0.74 * l(k,j) / ll(k,j) ) * &
	179	e(k,j,i) * SQRT( e(k,j,i) ) / l(k,j)
[1]	180
[65]	181	tend(k,j,i) = tend(k,j,i) &
[1]	182	+ ( &
	183	( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) ) &
	184	- ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) ) &
	185	) * ddx2 &
	186	+ ( &
	187	( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) ) &
	188	- ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) ) &
	189	) * ddy2 &
	190	+ ( &
	191	( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
	192	- ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k) &
	193	) * ddzw(k) &
	194	- dissipation(k,j)
	195
[65]	196	ENDDO
[1]	197	ENDDO
[65]	198
	199	!
	200	!-- Store dissipation if needed for calculating the sgs particle
	201	!-- velocities
[1065]	202	IF ( use_sgs_for_particles .OR. wang_kernel .OR. &
	203	turbulence ) THEN
[65]	204	DO j = nys, nyn
	205	DO k = nzb_s_inner(j,i)+1, nzt
	206	diss(k,j,i) = dissipation(k,j)
	207	ENDDO
	208	ENDDO
	209	ENDIF
	210
[1]	211	ENDDO
	212
[65]	213	ELSE
	214
	215	DO i = nxl, nxr
	216	DO j = nys, nyn
	217	DO k = nzb_s_inner(j,i)+1, nzt
	218	!
	219	!-- Calculate the mixing length (for dissipation)
[97]	220	dvar_dz = atmos_ocean_sign * &
	221	( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
	222	IF ( dvar_dz > 0.0 ) THEN
[65]	223	l_stable = 0.76 * SQRT( e(k,j,i) ) / &
[97]	224	SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5
[65]	225	ELSE
	226	l_stable = l_grid(k)
	227	ENDIF
	228	!
	229	!-- Adjustment of the mixing length
	230	IF ( wall_adjustment ) THEN
[94]	231	l(k,j) = MIN( wall_adjustment_factor * &
	232	( zu(k) - zw(nzb_s_inner(j,i)) ), &
	233	l_grid(k), l_stable )
	234	ll(k,j) = MIN( wall_adjustment_factor * &
	235	( zu(k) - zw(nzb_s_inner(j,i)) ), &
	236	l_grid(k) )
[65]	237	ELSE
	238	l(k,j) = MIN( l_grid(k), l_stable )
	239	ll(k,j) = l_grid(k)
	240	ENDIF
	241
	242	ENDDO
	243	ENDDO
	244
	245	!
	246	!-- Calculate the tendency terms
[1]	247	DO j = nys, nyn
[19]	248	DO k = nzb_s_inner(j,i)+1, nzt
[65]	249
	250	dissipation(k,j) = ( 0.19 + 0.74 * l(k,j) / ll(k,j) ) * &
	251	e(k,j,i) * SQRT( e(k,j,i) ) / l(k,j)
	252
	253	tend(k,j,i) = tend(k,j,i) &
	254	+ ( &
	255	( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) ) &
	256	- ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) ) &
	257	) * ddx2 &
	258	+ ( &
	259	( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) ) &
	260	- ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) ) &
	261	) * ddy2 &
	262	+ ( &
	263	( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
	264	- ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k) &
	265	) * ddzw(k) &
	266	- dissipation(k,j)
	267
[1]	268	ENDDO
	269	ENDDO
	270
[65]	271	!
	272	!-- Store dissipation if needed for calculating the sgs particle
	273	!-- velocities
[1065]	274	IF ( use_sgs_for_particles .OR. wang_kernel .OR. &
	275	turbulence ) THEN
[65]	276	DO j = nys, nyn
	277	DO k = nzb_s_inner(j,i)+1, nzt
	278	diss(k,j,i) = dissipation(k,j)
	279	ENDDO
	280	ENDDO
	281	ENDIF
[1]	282
[65]	283	ENDDO
	284
	285	ENDIF
	286
[1]	287	!
	288	!-- Boundary condition for dissipation
[1065]	289	IF ( use_sgs_for_particles .OR. wang_kernel .OR. turbulence ) THEN
[1]	290	DO i = nxl, nxr
	291	DO j = nys, nyn
	292	diss(nzb_s_inner(j,i),j,i) = diss(nzb_s_inner(j,i)+1,j,i)
	293	ENDDO
	294	ENDDO
	295	ENDIF
	296
	297	END SUBROUTINE diffusion_e
	298
	299
	300	!------------------------------------------------------------------------------!
[1015]	301	! Call for all grid points - accelerator version
	302	!------------------------------------------------------------------------------!
	303	SUBROUTINE diffusion_e_acc( var, var_reference )
	304
[1320]	305	USE arrays_3d, &
	306	ONLY: dd2zu, ddzu, ddzw, diss, e, km, l_grid, tend, zu, zw
	307
	308	USE control_parameters, &
	309	ONLY: atmos_ocean_sign, g, turbulence, use_single_reference_value, &
	310	wall_adjustment, wall_adjustment_factor
	311
	312	USE grid_variables, &
	313	ONLY: ddx2, ddy2
	314
	315	USE indices, &
	316	ONLY: i_left, i_right, j_north, j_south, nzb_s_inner, nzt
	317
	318	USE kinds
	319
	320	USE particle_attributes, &
	321	ONLY: use_sgs_for_particles, wang_kernel
[1015]	322
	323	IMPLICIT NONE
	324
[1320]	325	INTEGER(iwp) :: i !:
	326	INTEGER(iwp) :: j !:
	327	INTEGER(iwp) :: k !:
	328	REAL(wp) :: dissipation !:
	329	REAL(wp) :: dvar_dz !:
	330	REAL(wp) :: l !:
	331	REAL(wp) :: ll !:
	332	REAL(wp) :: l_stable !:
	333	REAL(wp) :: var_reference !:
[1015]	334
	335	#if defined( __nopointer )
[1320]	336	REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: var !:
[1015]	337	#else
[1320]	338	REAL(wp), DIMENSION(:,:,:), POINTER :: var !:
[1015]	339	#endif
	340
	341
	342	!
	343	!-- This if clause must be outside the k-loop because otherwise
	344	!-- runtime errors occur with -C hopt on NEC
[1179]	345	IF ( use_single_reference_value ) THEN
[1171]	346
	347	!$acc kernels present( ddzu, ddzw, dd2zu, diss, e, km, l_grid ) &
	348	!$acc present( nzb_s_inner, rif, tend, var, zu, zw )
	349	DO i = i_left, i_right
	350	DO j = j_south, j_north
	351	DO k = 1, nzt
	352
	353	IF ( k > nzb_s_inner(j,i) ) THEN
[1015]	354	!
[1171]	355	!-- Calculate the mixing length (for dissipation)
	356	dvar_dz = atmos_ocean_sign * &
	357	( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
	358	IF ( dvar_dz > 0.0 ) THEN
	359	l_stable = 0.76 * SQRT( e(k,j,i) ) / &
	360	SQRT( g / var_reference * dvar_dz ) + 1E-5
	361	ELSE
	362	l_stable = l_grid(k)
	363	ENDIF
[1015]	364	!
[1171]	365	!-- Adjustment of the mixing length
	366	IF ( wall_adjustment ) THEN
	367	l = MIN( wall_adjustment_factor * &
	368	( zu(k) - zw(nzb_s_inner(j,i)) ), &
	369	l_grid(k), l_stable )
	370	ll = MIN( wall_adjustment_factor * &
	371	( zu(k) - zw(nzb_s_inner(j,i)) ), &
	372	l_grid(k) )
	373	ELSE
	374	l = MIN( l_grid(k), l_stable )
	375	ll = l_grid(k)
	376	ENDIF
[1015]	377	!
[1171]	378	!-- Calculate the tendency terms
	379	dissipation = ( 0.19 + 0.74 * l / ll ) * &
	380	e(k,j,i) * SQRT( e(k,j,i) ) / l
	381
	382	tend(k,j,i) = tend(k,j,i) &
	383	+ ( &
	384	( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) ) &
	385	- ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) ) &
	386	) * ddx2 &
	387	+ ( &
	388	( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) ) &
	389	- ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) ) &
	390	) * ddy2 &
	391	+ ( &
	392	( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
	393	- ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k) &
	394	) * ddzw(k) &
	395	- dissipation
	396
[1015]	397	!
[1171]	398	!-- Store dissipation if needed for calculating the sgs particle
	399	!-- velocities
	400	IF ( use_sgs_for_particles .OR. wang_kernel .OR. &
	401	turbulence ) THEN
	402	diss(k,j,i) = dissipation
	403	ENDIF
	404
	405	ENDIF
	406
	407	ENDDO
	408	ENDDO
	409	ENDDO
	410	!$acc end kernels
	411
[1015]	412	ELSE
	413
	414	!$acc kernels present( ddzu, ddzw, dd2zu, diss, e, km, l_grid ) &
	415	!$acc present( nzb_s_inner, rif, tend, var, zu, zw )
[1128]	416	DO i = i_left, i_right
	417	DO j = j_south, j_north
[1015]	418	DO k = 1, nzt
	419
	420	IF ( k > nzb_s_inner(j,i) ) THEN
	421	!
	422	!-- Calculate the mixing length (for dissipation)
	423	dvar_dz = atmos_ocean_sign * &
	424	( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
	425	IF ( dvar_dz > 0.0 ) THEN
	426	l_stable = 0.76 * SQRT( e(k,j,i) ) / &
	427	SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5
	428	ELSE
	429	l_stable = l_grid(k)
	430	ENDIF
	431	!
	432	!-- Adjustment of the mixing length
	433	IF ( wall_adjustment ) THEN
	434	l = MIN( wall_adjustment_factor * &
	435	( zu(k) - zw(nzb_s_inner(j,i)) ), &
	436	l_grid(k), l_stable )
	437	ll = MIN( wall_adjustment_factor * &
	438	( zu(k) - zw(nzb_s_inner(j,i)) ), &
	439	l_grid(k) )
	440	ELSE
	441	l = MIN( l_grid(k), l_stable )
	442	ll = l_grid(k)
	443	ENDIF
	444	!
	445	!-- Calculate the tendency terms
	446	dissipation = ( 0.19 + 0.74 * l / ll ) * &
	447	e(k,j,i) * SQRT( e(k,j,i) ) / l
	448
	449	tend(k,j,i) = tend(k,j,i) &
	450	+ ( &
	451	( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) ) &
	452	- ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) ) &
	453	) * ddx2 &
	454	+ ( &
	455	( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) ) &
	456	- ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) ) &
	457	) * ddy2 &
	458	+ ( &
	459	( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
	460	- ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k) &
	461	) * ddzw(k) &
[1171]	462	- dissipation
[1015]	463
	464	!
	465	!-- Store dissipation if needed for calculating the sgs
	466	!-- particle velocities
[1065]	467	IF ( use_sgs_for_particles .OR. wang_kernel .OR. &
	468	turbulence ) THEN
[1015]	469	diss(k,j,i) = dissipation
	470	ENDIF
	471
	472	ENDIF
	473
	474	ENDDO
	475	ENDDO
	476	ENDDO
	477	!$acc end kernels
	478
	479	ENDIF
	480
	481	!
	482	!-- Boundary condition for dissipation
[1065]	483	IF ( use_sgs_for_particles .OR. wang_kernel .OR. turbulence ) THEN
[1015]	484	!$acc kernels present( diss, nzb_s_inner )
[1128]	485	DO i = i_left, i_right
	486	DO j = j_south, j_north
[1015]	487	diss(nzb_s_inner(j,i),j,i) = diss(nzb_s_inner(j,i)+1,j,i)
	488	ENDDO
	489	ENDDO
	490	!$acc end kernels
	491	ENDIF
	492
	493	END SUBROUTINE diffusion_e_acc
	494
	495
	496	!------------------------------------------------------------------------------!
[1]	497	! Call for grid point i,j
	498	!------------------------------------------------------------------------------!
[1001]	499	SUBROUTINE diffusion_e_ij( i, j, var, var_reference )
[1]	500
[1320]	501	USE arrays_3d, &
	502	ONLY: dd2zu, ddzu, ddzw, diss, e, km, l_grid, tend, zu, zw
	503
	504	USE control_parameters, &
	505	ONLY: atmos_ocean_sign, g, turbulence, use_single_reference_value, &
	506	wall_adjustment, wall_adjustment_factor
	507
	508	USE grid_variables, &
	509	ONLY: ddx2, ddy2
	510
	511	USE indices, &
	512	ONLY: nzb, nzb_s_inner, nzt
	513
	514	USE kinds
	515
	516	USE particle_attributes, &
	517	ONLY: use_sgs_for_particles, wang_kernel
[1]	518
	519	IMPLICIT NONE
	520
[1320]	521	INTEGER(iwp) :: i !:
	522	INTEGER(iwp) :: j !:
	523	INTEGER(iwp) :: k !:
	524	REAL(wp) :: dvar_dz !:
	525	REAL(wp) :: l_stable !:
	526	REAL(wp) :: var_reference !:
[1]	527
[1010]	528	#if defined( __nopointer )
[1320]	529	REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: var !:
[1010]	530	#else
[1320]	531	REAL(wp), DIMENSION(:,:,:), POINTER :: var !:
[1010]	532	#endif
[1320]	533	REAL(wp), DIMENSION(nzb+1:nzt) :: dissipation !:
	534	REAL(wp), DIMENSION(nzb+1:nzt) :: l !:
	535	REAL(wp), DIMENSION(nzb+1:nzt) :: ll !:
[1]	536
[1001]	537
[1]	538	!
	539	!-- Calculate the mixing length (for dissipation)
[19]	540	DO k = nzb_s_inner(j,i)+1, nzt
[97]	541	dvar_dz = atmos_ocean_sign * &
	542	( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
	543	IF ( dvar_dz > 0.0 ) THEN
[1179]	544	IF ( use_single_reference_value ) THEN
[57]	545	l_stable = 0.76 * SQRT( e(k,j,i) ) / &
[97]	546	SQRT( g / var_reference * dvar_dz ) + 1E-5
[57]	547	ELSE
	548	l_stable = 0.76 * SQRT( e(k,j,i) ) / &
[97]	549	SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5
[57]	550	ENDIF
[1]	551	ELSE
	552	l_stable = l_grid(k)
	553	ENDIF
	554	!
	555	!-- Adjustment of the mixing length
	556	IF ( wall_adjustment ) THEN
[94]	557	l(k) = MIN( wall_adjustment_factor * &
	558	( zu(k) - zw(nzb_s_inner(j,i)) ), l_grid(k), &
	559	l_stable )
	560	ll(k) = MIN( wall_adjustment_factor * &
	561	( zu(k) - zw(nzb_s_inner(j,i)) ), l_grid(k) )
[1]	562	ELSE
	563	l(k) = MIN( l_grid(k), l_stable )
	564	ll(k) = l_grid(k)
	565	ENDIF
	566	!
	567	!-- Calculate the tendency term
	568	dissipation(k) = ( 0.19 + 0.74 * l(k) / ll(k) ) * e(k,j,i) * &
	569	SQRT( e(k,j,i) ) / l(k)
	570
	571	tend(k,j,i) = tend(k,j,i) &
	572	+ ( &
	573	( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) ) &
	574	- ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) ) &
	575	) * ddx2 &
	576	+ ( &
	577	( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) ) &
	578	- ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) ) &
	579	) * ddy2 &
	580	+ ( &
	581	( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
	582	- ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k) &
	583	) * ddzw(k) &
	584	- dissipation(k)
	585
	586	ENDDO
	587
	588	!
	589	!-- Store dissipation if needed for calculating the sgs particle velocities
[1065]	590	IF ( use_sgs_for_particles .OR. wang_kernel .OR. turbulence ) THEN
[19]	591	DO k = nzb_s_inner(j,i)+1, nzt
[1]	592	diss(k,j,i) = dissipation(k)
	593	ENDDO
	594	!
	595	!-- Boundary condition for dissipation
	596	diss(nzb_s_inner(j,i),j,i) = diss(nzb_s_inner(j,i)+1,j,i)
	597	ENDIF
	598
	599	END SUBROUTINE diffusion_e_ij
	600
	601	END MODULE diffusion_e_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |