Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |