Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |