Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |