Home

Context Navigation

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

Last change on this file since 1320 was 1320, checked in by raasch, 10 years ago
ONLY-attribute added to USE-statements, kind-parameters added to all INTEGER and REAL declaration statements, kinds are defined in new module kinds, old module precision_kind is removed, revision history before 2012 removed, comment fields (!:) to be used for variable explanations added to all variable declaration statements
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	! -----------------
[1320]	22	! ONLY-attribute added to USE-statements,
	23	! kind-parameters added to all INTEGER and REAL declaration statements,
	24	! kinds are defined in new module kinds,
	25	! old module precision_kind is removed,
	26	! revision history before 2012 removed,
	27	! comment fields (!:) to be used for variable explanations added to
	28	! all variable declaration statements
[98]	29	!
	30	! Former revisions:
	31	! -----------------
	32	! $Id: diffusion_e.f90 1320 2014-03-20 08:40:49Z raasch $
	33	!
[1258]	34	! 1257 2013-11-08 15:18:40Z raasch
	35	! openacc loop and loop vector clauses removed
	36	!
[1182]	37	! 1179 2013-06-14 05:57:58Z raasch
	38	! use_reference renamed use_single_reference_value
	39	!
[1172]	40	! 1171 2013-05-30 11:27:45Z raasch
	41	! use_reference-case activated in accelerator version
	42	!
[1132]	43	! 1128 2013-04-12 06:19:32Z raasch
	44	! loop index bounds in accelerator version replaced by i_left, i_right, j_south,
	45	! j_north
	46	!
[1066]	47	! 1065 2012-11-22 17:42:36Z hoffmann
	48	! Enabled the claculation of diss in case of turbulence = .TRUE. (parameterized
	49	! effects of turbulence on autoconversion and accretion in two-moments cloud
	50	! physics scheme).
	51	!
[1037]	52	! 1036 2012-10-22 13:43:42Z raasch
	53	! code put under GPL (PALM 3.9)
	54	!
[1017]	55	! 1015 2012-09-27 09:23:24Z raasch
	56	! accelerator version (*_acc) added,
	57	! adjustment of mixing length to the Prandtl mixing length at first grid point
	58	! above ground removed
	59	!
[1011]	60	! 1010 2012-09-20 07:59:54Z raasch
	61	! cpp switch __nopointer added for pointer free version
	62	!
[1002]	63	! 1001 2012-09-13 14:08:46Z raasch
	64	! most arrays comunicated by module instead of parameter list
	65	!
[826]	66	! 825 2012-02-19 03:03:44Z raasch
	67	! wang_collision_kernel renamed wang_kernel
	68	!
[1]	69	! Revision 1.1 1997/09/19 07:40:24 raasch
	70	! Initial revision
	71	!
	72	!
	73	! Description:
	74	! ------------
	75	! Diffusion- and dissipation terms for the TKE
	76	!------------------------------------------------------------------------------!
	77
	78	PRIVATE
[1015]	79	PUBLIC diffusion_e, diffusion_e_acc
[1]	80
	81
	82	INTERFACE diffusion_e
	83	MODULE PROCEDURE diffusion_e
	84	MODULE PROCEDURE diffusion_e_ij
	85	END INTERFACE diffusion_e
	86
[1015]	87	INTERFACE diffusion_e_acc
	88	MODULE PROCEDURE diffusion_e_acc
	89	END INTERFACE diffusion_e_acc
	90
[1]	91	CONTAINS
	92
	93
	94	!------------------------------------------------------------------------------!
	95	! Call for all grid points
	96	!------------------------------------------------------------------------------!
[1001]	97	SUBROUTINE diffusion_e( var, var_reference )
[1]	98
[1320]	99	USE arrays_3d, &
	100	ONLY: dd2zu, ddzu, ddzw, diss, e, km, l_grid, tend, zu, zw
	101
	102	USE control_parameters, &
	103	ONLY: atmos_ocean_sign, g, turbulence, use_single_reference_value, &
	104	wall_adjustment, wall_adjustment_factor
	105
	106	USE grid_variables, &
	107	ONLY: ddx2, ddy2
	108
	109	USE indices, &
	110	ONLY: nxl, nxr, nyn, nys, nzb, nzb_s_inner, nzt
	111
	112	USE kinds
	113
	114	USE particle_attributes, &
	115	ONLY: use_sgs_for_particles, wang_kernel
[1]	116
	117	IMPLICIT NONE
	118
[1320]	119	INTEGER(iwp) :: i !:
	120	INTEGER(iwp) :: j !:
	121	INTEGER(iwp) :: k !:
	122	REAL(wp) :: dvar_dz !:
	123	REAL(wp) :: l_stable !:
	124	REAL(wp) :: var_reference !:
[1001]	125
[1010]	126	#if defined( __nopointer )
[1320]	127	REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: var !:
[1010]	128	#else
[1320]	129	REAL(wp), DIMENSION(:,:,:), POINTER :: var !:
[1010]	130	#endif
[1320]	131	REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: dissipation !:
	132	REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: l !:
	133	REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: 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
[1320]	303	USE arrays_3d, &
	304	ONLY: dd2zu, ddzu, ddzw, diss, e, km, l_grid, tend, zu, zw
	305
	306	USE control_parameters, &
	307	ONLY: atmos_ocean_sign, g, turbulence, use_single_reference_value, &
	308	wall_adjustment, wall_adjustment_factor
	309
	310	USE grid_variables, &
	311	ONLY: ddx2, ddy2
	312
	313	USE indices, &
	314	ONLY: i_left, i_right, j_north, j_south, nzb_s_inner, nzt
	315
	316	USE kinds
	317
	318	USE particle_attributes, &
	319	ONLY: use_sgs_for_particles, wang_kernel
[1015]	320
	321	IMPLICIT NONE
	322
[1320]	323	INTEGER(iwp) :: i !:
	324	INTEGER(iwp) :: j !:
	325	INTEGER(iwp) :: k !:
	326	REAL(wp) :: dissipation !:
	327	REAL(wp) :: dvar_dz !:
	328	REAL(wp) :: l !:
	329	REAL(wp) :: ll !:
	330	REAL(wp) :: l_stable !:
	331	REAL(wp) :: var_reference !:
[1015]	332
	333	#if defined( __nopointer )
[1320]	334	REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: var !:
[1015]	335	#else
[1320]	336	REAL(wp), DIMENSION(:,:,:), POINTER :: var !:
[1015]	337	#endif
	338
	339
	340	!
	341	!-- This if clause must be outside the k-loop because otherwise
	342	!-- runtime errors occur with -C hopt on NEC
[1179]	343	IF ( use_single_reference_value ) THEN
[1171]	344
	345	!$acc kernels present( ddzu, ddzw, dd2zu, diss, e, km, l_grid ) &
	346	!$acc present( nzb_s_inner, rif, tend, var, zu, zw )
	347	DO i = i_left, i_right
	348	DO j = j_south, j_north
	349	DO k = 1, nzt
	350
	351	IF ( k > nzb_s_inner(j,i) ) THEN
[1015]	352	!
[1171]	353	!-- Calculate the mixing length (for dissipation)
	354	dvar_dz = atmos_ocean_sign * &
	355	( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
	356	IF ( dvar_dz > 0.0 ) THEN
	357	l_stable = 0.76 * SQRT( e(k,j,i) ) / &
	358	SQRT( g / var_reference * dvar_dz ) + 1E-5
	359	ELSE
	360	l_stable = l_grid(k)
	361	ENDIF
[1015]	362	!
[1171]	363	!-- Adjustment of the mixing length
	364	IF ( wall_adjustment ) THEN
	365	l = MIN( wall_adjustment_factor * &
	366	( zu(k) - zw(nzb_s_inner(j,i)) ), &
	367	l_grid(k), l_stable )
	368	ll = MIN( wall_adjustment_factor * &
	369	( zu(k) - zw(nzb_s_inner(j,i)) ), &
	370	l_grid(k) )
	371	ELSE
	372	l = MIN( l_grid(k), l_stable )
	373	ll = l_grid(k)
	374	ENDIF
[1015]	375	!
[1171]	376	!-- Calculate the tendency terms
	377	dissipation = ( 0.19 + 0.74 * l / ll ) * &
	378	e(k,j,i) * SQRT( e(k,j,i) ) / l
	379
	380	tend(k,j,i) = tend(k,j,i) &
	381	+ ( &
	382	( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) ) &
	383	- ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) ) &
	384	) * ddx2 &
	385	+ ( &
	386	( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) ) &
	387	- ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) ) &
	388	) * ddy2 &
	389	+ ( &
	390	( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
	391	- ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k) &
	392	) * ddzw(k) &
	393	- dissipation
	394
[1015]	395	!
[1171]	396	!-- Store dissipation if needed for calculating the sgs particle
	397	!-- velocities
	398	IF ( use_sgs_for_particles .OR. wang_kernel .OR. &
	399	turbulence ) THEN
	400	diss(k,j,i) = dissipation
	401	ENDIF
	402
	403	ENDIF
	404
	405	ENDDO
	406	ENDDO
	407	ENDDO
	408	!$acc end kernels
	409
[1015]	410	ELSE
	411
	412	!$acc kernels present( ddzu, ddzw, dd2zu, diss, e, km, l_grid ) &
	413	!$acc present( nzb_s_inner, rif, tend, var, zu, zw )
[1128]	414	DO i = i_left, i_right
	415	DO j = j_south, j_north
[1015]	416	DO k = 1, nzt
	417
	418	IF ( k > nzb_s_inner(j,i) ) THEN
	419	!
	420	!-- Calculate the mixing length (for dissipation)
	421	dvar_dz = atmos_ocean_sign * &
	422	( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
	423	IF ( dvar_dz > 0.0 ) THEN
	424	l_stable = 0.76 * SQRT( e(k,j,i) ) / &
	425	SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5
	426	ELSE
	427	l_stable = l_grid(k)
	428	ENDIF
	429	!
	430	!-- Adjustment of the mixing length
	431	IF ( wall_adjustment ) THEN
	432	l = MIN( wall_adjustment_factor * &
	433	( zu(k) - zw(nzb_s_inner(j,i)) ), &
	434	l_grid(k), l_stable )
	435	ll = MIN( wall_adjustment_factor * &
	436	( zu(k) - zw(nzb_s_inner(j,i)) ), &
	437	l_grid(k) )
	438	ELSE
	439	l = MIN( l_grid(k), l_stable )
	440	ll = l_grid(k)
	441	ENDIF
	442	!
	443	!-- Calculate the tendency terms
	444	dissipation = ( 0.19 + 0.74 * l / ll ) * &
	445	e(k,j,i) * SQRT( e(k,j,i) ) / l
	446
	447	tend(k,j,i) = tend(k,j,i) &
	448	+ ( &
	449	( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) ) &
	450	- ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) ) &
	451	) * ddx2 &
	452	+ ( &
	453	( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) ) &
	454	- ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) ) &
	455	) * ddy2 &
	456	+ ( &
	457	( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
	458	- ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k) &
	459	) * ddzw(k) &
[1171]	460	- dissipation
[1015]	461
	462	!
	463	!-- Store dissipation if needed for calculating the sgs
	464	!-- particle velocities
[1065]	465	IF ( use_sgs_for_particles .OR. wang_kernel .OR. &
	466	turbulence ) THEN
[1015]	467	diss(k,j,i) = dissipation
	468	ENDIF
	469
	470	ENDIF
	471
	472	ENDDO
	473	ENDDO
	474	ENDDO
	475	!$acc end kernels
	476
	477	ENDIF
	478
	479	!
	480	!-- Boundary condition for dissipation
[1065]	481	IF ( use_sgs_for_particles .OR. wang_kernel .OR. turbulence ) THEN
[1015]	482	!$acc kernels present( diss, nzb_s_inner )
[1128]	483	DO i = i_left, i_right
	484	DO j = j_south, j_north
[1015]	485	diss(nzb_s_inner(j,i),j,i) = diss(nzb_s_inner(j,i)+1,j,i)
	486	ENDDO
	487	ENDDO
	488	!$acc end kernels
	489	ENDIF
	490
	491	END SUBROUTINE diffusion_e_acc
	492
	493
	494	!------------------------------------------------------------------------------!
[1]	495	! Call for grid point i,j
	496	!------------------------------------------------------------------------------!
[1001]	497	SUBROUTINE diffusion_e_ij( i, j, var, var_reference )
[1]	498
[1320]	499	USE arrays_3d, &
	500	ONLY: dd2zu, ddzu, ddzw, diss, e, km, l_grid, tend, zu, zw
	501
	502	USE control_parameters, &
	503	ONLY: atmos_ocean_sign, g, turbulence, use_single_reference_value, &
	504	wall_adjustment, wall_adjustment_factor
	505
	506	USE grid_variables, &
	507	ONLY: ddx2, ddy2
	508
	509	USE indices, &
	510	ONLY: nzb, nzb_s_inner, nzt
	511
	512	USE kinds
	513
	514	USE particle_attributes, &
	515	ONLY: use_sgs_for_particles, wang_kernel
[1]	516
	517	IMPLICIT NONE
	518
[1320]	519	INTEGER(iwp) :: i !:
	520	INTEGER(iwp) :: j !:
	521	INTEGER(iwp) :: k !:
	522	REAL(wp) :: dvar_dz !:
	523	REAL(wp) :: l_stable !:
	524	REAL(wp) :: var_reference !:
[1]	525
[1010]	526	#if defined( __nopointer )
[1320]	527	REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: var !:
[1010]	528	#else
[1320]	529	REAL(wp), DIMENSION(:,:,:), POINTER :: var !:
[1010]	530	#endif
[1320]	531	REAL(wp), DIMENSION(nzb+1:nzt) :: dissipation !:
	532	REAL(wp), DIMENSION(nzb+1:nzt) :: l !:
	533	REAL(wp), DIMENSION(nzb+1:nzt) :: ll !:
[1]	534
[1001]	535
[1]	536	!
	537	!-- Calculate the mixing length (for dissipation)
[19]	538	DO k = nzb_s_inner(j,i)+1, nzt
[97]	539	dvar_dz = atmos_ocean_sign * &
	540	( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
	541	IF ( dvar_dz > 0.0 ) THEN
[1179]	542	IF ( use_single_reference_value ) THEN
[57]	543	l_stable = 0.76 * SQRT( e(k,j,i) ) / &
[97]	544	SQRT( g / var_reference * dvar_dz ) + 1E-5
[57]	545	ELSE
	546	l_stable = 0.76 * SQRT( e(k,j,i) ) / &
[97]	547	SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5
[57]	548	ENDIF
[1]	549	ELSE
	550	l_stable = l_grid(k)
	551	ENDIF
	552	!
	553	!-- Adjustment of the mixing length
	554	IF ( wall_adjustment ) THEN
[94]	555	l(k) = MIN( wall_adjustment_factor * &
	556	( zu(k) - zw(nzb_s_inner(j,i)) ), l_grid(k), &
	557	l_stable )
	558	ll(k) = MIN( wall_adjustment_factor * &
	559	( zu(k) - zw(nzb_s_inner(j,i)) ), l_grid(k) )
[1]	560	ELSE
	561	l(k) = MIN( l_grid(k), l_stable )
	562	ll(k) = l_grid(k)
	563	ENDIF
	564	!
	565	!-- Calculate the tendency term
	566	dissipation(k) = ( 0.19 + 0.74 * l(k) / ll(k) ) * e(k,j,i) * &
	567	SQRT( e(k,j,i) ) / l(k)
	568
	569	tend(k,j,i) = tend(k,j,i) &
	570	+ ( &
	571	( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) ) &
	572	- ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) ) &
	573	) * ddx2 &
	574	+ ( &
	575	( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) ) &
	576	- ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) ) &
	577	) * ddy2 &
	578	+ ( &
	579	( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
	580	- ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k) &
	581	) * ddzw(k) &
	582	- dissipation(k)
	583
	584	ENDDO
	585
	586	!
	587	!-- Store dissipation if needed for calculating the sgs particle velocities
[1065]	588	IF ( use_sgs_for_particles .OR. wang_kernel .OR. turbulence ) THEN
[19]	589	DO k = nzb_s_inner(j,i)+1, nzt
[1]	590	diss(k,j,i) = dissipation(k)
	591	ENDDO
	592	!
	593	!-- Boundary condition for dissipation
	594	diss(nzb_s_inner(j,i),j,i) = diss(nzb_s_inner(j,i)+1,j,i)
	595	ENDIF
	596
	597	END SUBROUTINE diffusion_e_ij
	598
	599	END MODULE diffusion_e_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |