Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |