Home

Context Navigation

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

Last change on this file since 2037 was 2037, checked in by knoop, 7 years ago
Anelastic approximation implemented
Property svn:keywords set to `Id`
File size: 26.5 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |