Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |