Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |