Home

Context Navigation

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

Last change on this file since 1806 was 1683, checked in by knoop, 9 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 24.7 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |