Home

Context Navigation

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

Last change on this file since 2412 was 2233, checked in by suehring, 7 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 19.7 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	!
[2101]	17	! Copyright 1997-2017 Leibniz Universitaet Hannover
[2000]	18	!------------------------------------------------------------------------------!
[1036]	19	!
[484]	20	! Current revisions:
[1]	21	! -----------------
[1341]	22	!
[2233]	23	!
[1321]	24	! Former revisions:
	25	! -----------------
	26	! $Id: diffusion_e.f90 2233 2017-05-30 18:08:54Z raasch $
	27	!
[2233]	28	! 2232 2017-05-30 17:47:52Z suehring
	29	! Adjustments to new topography and surface concept
	30	!
[2119]	31	! 2118 2017-01-17 16:38:49Z raasch
	32	! OpenACC version of subroutine removed
	33	!
[2038]	34	! 2037 2016-10-26 11:15:40Z knoop
	35	! Anelastic approximation implemented
	36	!
[2001]	37	! 2000 2016-08-20 18:09:15Z knoop
	38	! Forced header and separation lines into 80 columns
	39	!
[1874]	40	! 1873 2016-04-18 14:50:06Z maronga
	41	! Module renamed (removed _mod)
	42	!
[1851]	43	! 1850 2016-04-08 13:29:27Z maronga
	44	! Module renamed
	45	! Adapted for modularization of microphysics
	46	!
[1832]	47	! 1831 2016-04-07 13:15:51Z hoffmann
	48	! turbulence renamed collision_turbulence
	49	!
[1683]	50	! 1682 2015-10-07 23:56:08Z knoop
	51	! Code annotations made doxygen readable
	52	!
[1375]	53	! 1374 2014-04-25 12:55:07Z raasch
	54	! missing variables added to ONLY list
	55	! rif removed from acc-present-list
	56	!
[1341]	57	! 1340 2014-03-25 19:45:13Z kanani
	58	! REAL constants defined as wp-kind
	59	!
[1321]	60	! 1320 2014-03-20 08:40:49Z raasch
[1320]	61	! ONLY-attribute added to USE-statements,
	62	! kind-parameters added to all INTEGER and REAL declaration statements,
	63	! kinds are defined in new module kinds,
	64	! revision history before 2012 removed,
	65	! comment fields (!:) to be used for variable explanations added to
	66	! all variable declaration statements
[98]	67	!
[1258]	68	! 1257 2013-11-08 15:18:40Z raasch
	69	! openacc loop and loop vector clauses removed
	70	!
[1182]	71	! 1179 2013-06-14 05:57:58Z raasch
	72	! use_reference renamed use_single_reference_value
	73	!
[1172]	74	! 1171 2013-05-30 11:27:45Z raasch
	75	! use_reference-case activated in accelerator version
	76	!
[1132]	77	! 1128 2013-04-12 06:19:32Z raasch
	78	! loop index bounds in accelerator version replaced by i_left, i_right, j_south,
	79	! j_north
	80	!
[1066]	81	! 1065 2012-11-22 17:42:36Z hoffmann
	82	! Enabled the claculation of diss in case of turbulence = .TRUE. (parameterized
	83	! effects of turbulence on autoconversion and accretion in two-moments cloud
	84	! physics scheme).
	85	!
[1037]	86	! 1036 2012-10-22 13:43:42Z raasch
	87	! code put under GPL (PALM 3.9)
	88	!
[1017]	89	! 1015 2012-09-27 09:23:24Z raasch
	90	! accelerator version (*_acc) added,
	91	! adjustment of mixing length to the Prandtl mixing length at first grid point
	92	! above ground removed
	93	!
[1011]	94	! 1010 2012-09-20 07:59:54Z raasch
	95	! cpp switch __nopointer added for pointer free version
	96	!
[1002]	97	! 1001 2012-09-13 14:08:46Z raasch
	98	! most arrays comunicated by module instead of parameter list
	99	!
[826]	100	! 825 2012-02-19 03:03:44Z raasch
	101	! wang_collision_kernel renamed wang_kernel
	102	!
[1]	103	! Revision 1.1 1997/09/19 07:40:24 raasch
	104	! Initial revision
	105	!
	106	!
	107	! Description:
	108	! ------------
[1682]	109	!> Diffusion- and dissipation terms for the TKE
[1]	110	!------------------------------------------------------------------------------!
[1682]	111	MODULE diffusion_e_mod
	112
[1]	113
	114	PRIVATE
[2118]	115	PUBLIC diffusion_e
[1]	116
	117
	118	INTERFACE diffusion_e
	119	MODULE PROCEDURE diffusion_e
	120	MODULE PROCEDURE diffusion_e_ij
	121	END INTERFACE diffusion_e
	122
	123	CONTAINS
	124
	125
	126	!------------------------------------------------------------------------------!
[1682]	127	! Description:
	128	! ------------
	129	!> Call for all grid points
[1]	130	!------------------------------------------------------------------------------!
[1001]	131	SUBROUTINE diffusion_e( var, var_reference )
[1]	132
[1320]	133	USE arrays_3d, &
[2232]	134	ONLY: dd2zu, ddzu, ddzw, diss, e, km, l_grid, l_wall, tend, zu, zw,&
[2037]	135	drho_air, rho_air_zw
[1320]	136
	137	USE control_parameters, &
[1831]	138	ONLY: atmos_ocean_sign, g, use_single_reference_value, &
[1320]	139	wall_adjustment, wall_adjustment_factor
[1831]	140
[1320]	141	USE grid_variables, &
	142	ONLY: ddx2, ddy2
	143
	144	USE indices, &
[2232]	145	ONLY: nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzb_max, &
	146	nzt, wall_flags_0
[1320]	147
	148	USE kinds
[1849]	149
	150	USE microphysics_mod, &
	151	ONLY: collision_turbulence
	152
[1320]	153	USE particle_attributes, &
	154	ONLY: use_sgs_for_particles, wang_kernel
[1]	155
[2232]	156	USE surface_mod, &
	157	ONLY : bc_h
	158
[1]	159	IMPLICIT NONE
	160
[2232]	161	INTEGER(iwp) :: i !< running index x direction
	162	INTEGER(iwp) :: j !< running index y direction
	163	INTEGER(iwp) :: k !< running index z direction
	164	INTEGER(iwp) :: m !< running index surface elements
	165
[1682]	166	REAL(wp) :: dvar_dz !<
[2232]	167	REAL(wp) :: flag !< flag to mask topography
[1682]	168	REAL(wp) :: l_stable !<
	169	REAL(wp) :: var_reference !<
[1001]	170
[1010]	171	#if defined( __nopointer )
[1682]	172	REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: var !<
[1010]	173	#else
[1682]	174	REAL(wp), DIMENSION(:,:,:), POINTER :: var !<
[1010]	175	#endif
[1682]	176	REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: dissipation !<
	177	REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: l !<
	178	REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: ll !<
[1]	179
	180
	181	!
[65]	182	!-- This if clause must be outside the k-loop because otherwise
	183	!-- runtime errors occur with -C hopt on NEC
[1179]	184	IF ( use_single_reference_value ) THEN
[65]	185
	186	DO i = nxl, nxr
	187	DO j = nys, nyn
[2232]	188	DO k = nzb+1, nzt
[1]	189	!
[65]	190	!-- Calculate the mixing length (for dissipation)
[97]	191	dvar_dz = atmos_ocean_sign * &
	192	( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
[1340]	193	IF ( dvar_dz > 0.0_wp ) THEN
	194	l_stable = 0.76_wp * SQRT( e(k,j,i) ) / &
	195	SQRT( g / var_reference * dvar_dz ) + 1E-5_wp
[57]	196	ELSE
[65]	197	l_stable = l_grid(k)
[57]	198	ENDIF
[1]	199	!
[65]	200	!-- Adjustment of the mixing length
	201	IF ( wall_adjustment ) THEN
[2232]	202	l(k,j) = MIN( wall_adjustment_factor * l_wall(k,j,i), &
[94]	203	l_grid(k), l_stable )
[2232]	204	ll(k,j) = MIN( wall_adjustment_factor * l_wall(k,j,i), &
[94]	205	l_grid(k) )
[65]	206	ELSE
	207	l(k,j) = MIN( l_grid(k), l_stable )
	208	ll(k,j) = l_grid(k)
	209	ENDIF
[1]	210
[65]	211	ENDDO
[1]	212	ENDDO
[65]	213
[1]	214	!
[65]	215	!-- Calculate the tendency terms
	216	DO j = nys, nyn
[2232]	217	DO k = nzb+1, nzt
	218	!
	219	!-- Predetermine flag to mask topography
	220	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
[1]	221
[2232]	222	dissipation(k,j) = ( 0.19_wp + 0.74_wp * l(k,j) / ll(k,j) ) * &
[65]	223	e(k,j,i) * SQRT( e(k,j,i) ) / l(k,j)
[1]	224
[2232]	225	tend(k,j,i) = tend(k,j,i) &
[1]	226	+ ( &
	227	( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) ) &
	228	- ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) ) &
[2232]	229	) * ddx2 * flag &
[1]	230	+ ( &
	231	( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) ) &
	232	- ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) ) &
[2232]	233	) * ddy2 * flag &
[1]	234	+ ( &
	235	( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
[2037]	236	* rho_air_zw(k) &
[1]	237	- ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k) &
[2037]	238	* rho_air_zw(k-1) &
[2232]	239	) * ddzw(k) * drho_air(k) * flag &
	240	- dissipation(k,j) * flag
[1]	241
[65]	242	ENDDO
[1]	243	ENDDO
[65]	244
	245	!
	246	!-- Store dissipation if needed for calculating the sgs particle
	247	!-- velocities
[1065]	248	IF ( use_sgs_for_particles .OR. wang_kernel .OR. &
[1831]	249	collision_turbulence ) THEN
[65]	250	DO j = nys, nyn
[2232]	251	DO k = nzb+1, nzt
	252	diss(k,j,i) = dissipation(k,j) * &
	253	MERGE( 1.0_wp, 0.0_wp, &
	254	BTEST( wall_flags_0(k,j,i), 0 ) )
[65]	255	ENDDO
	256	ENDDO
	257	ENDIF
	258
[1]	259	ENDDO
	260
[65]	261	ELSE
	262
	263	DO i = nxl, nxr
	264	DO j = nys, nyn
[2232]	265	DO k = nzb+1, nzt
[65]	266	!
	267	!-- Calculate the mixing length (for dissipation)
[97]	268	dvar_dz = atmos_ocean_sign * &
	269	( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
[1340]	270	IF ( dvar_dz > 0.0_wp ) THEN
	271	l_stable = 0.76_wp * SQRT( e(k,j,i) ) / &
	272	SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5_wp
[65]	273	ELSE
	274	l_stable = l_grid(k)
	275	ENDIF
	276	!
	277	!-- Adjustment of the mixing length
	278	IF ( wall_adjustment ) THEN
[2232]	279	l(k,j) = MIN( wall_adjustment_factor * l_wall(k,j,i), &
[94]	280	l_grid(k), l_stable )
[2232]	281	ll(k,j) = MIN( wall_adjustment_factor * l_wall(k,j,i), &
[94]	282	l_grid(k) )
[65]	283	ELSE
	284	l(k,j) = MIN( l_grid(k), l_stable )
	285	ll(k,j) = l_grid(k)
	286	ENDIF
	287
	288	ENDDO
	289	ENDDO
	290
	291	!
	292	!-- Calculate the tendency terms
[1]	293	DO j = nys, nyn
[2232]	294	DO k = nzb+1, nzt
	295	!
	296	!-- Predetermine flag to mask topography
	297	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
[65]	298
[2232]	299	dissipation(k,j) = ( 0.19_wp + 0.74_wp * l(k,j) / ll(k,j) ) *&
[1340]	300	e(k,j,i) * SQRT( e(k,j,i) ) / l(k,j)
[65]	301
[2232]	302	tend(k,j,i) = tend(k,j,i) &
[65]	303	+ ( &
	304	( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) ) &
	305	- ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) ) &
[2232]	306	) * ddx2 * flag &
[65]	307	+ ( &
	308	( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) ) &
	309	- ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) ) &
[2232]	310	) * ddy2 * flag &
[65]	311	+ ( &
	312	( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
[2037]	313	* rho_air_zw(k) &
[65]	314	- ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k) &
[2037]	315	* rho_air_zw(k-1) &
[2232]	316	) * ddzw(k) * drho_air(k) * flag &
	317	- dissipation(k,j) * flag
[65]	318
[1]	319	ENDDO
	320	ENDDO
	321
[65]	322	!
	323	!-- Store dissipation if needed for calculating the sgs particle
	324	!-- velocities
[1065]	325	IF ( use_sgs_for_particles .OR. wang_kernel .OR. &
[1831]	326	collision_turbulence ) THEN
[65]	327	DO j = nys, nyn
[2232]	328	DO k = nzb+1, nzt
	329	diss(k,j,i) = dissipation(k,j) * &
	330	MERGE( 1.0_wp, 0.0_wp, &
	331	BTEST( wall_flags_0(k,j,i), 0 ) )
[65]	332	ENDDO
	333	ENDDO
	334	ENDIF
[1]	335
[65]	336	ENDDO
	337
	338	ENDIF
	339
[1]	340	!
[2232]	341	!-- Neumann boundary condition for dissipation diss(nzb,:,:) = diss(nzb+1,:,:)
[1849]	342	IF ( use_sgs_for_particles .OR. wang_kernel .OR. &
[1831]	343	collision_turbulence ) THEN
[2232]	344	!
	345	!-- Upward facing surfaces
	346	DO m = 1, bc_h(0)%ns
	347	i = bc_h(0)%i(m)
	348	j = bc_h(0)%j(m)
	349	k = bc_h(0)%k(m)
	350	diss(k-1,j,i) = diss(k,j,i)
[1]	351	ENDDO
[2232]	352	!
	353	!-- Downward facing surfaces
	354	DO m = 1, bc_h(1)%ns
	355	i = bc_h(1)%i(m)
	356	j = bc_h(1)%j(m)
	357	k = bc_h(1)%k(m)
	358	diss(k+1,j,i) = diss(k,j,i)
	359	ENDDO
	360
[1]	361	ENDIF
	362
	363	END SUBROUTINE diffusion_e
	364
	365
	366	!------------------------------------------------------------------------------!
[1682]	367	! Description:
	368	! ------------
	369	!> Call for grid point i,j
[1]	370	!------------------------------------------------------------------------------!
[1001]	371	SUBROUTINE diffusion_e_ij( i, j, var, var_reference )
[1]	372
[1320]	373	USE arrays_3d, &
[2232]	374	ONLY: dd2zu, ddzu, ddzw, diss, e, km, l_grid, l_wall, tend, zu, zw,&
[2037]	375	drho_air, rho_air_zw
[1320]	376
	377	USE control_parameters, &
[1831]	378	ONLY: atmos_ocean_sign, g, use_single_reference_value, &
[1320]	379	wall_adjustment, wall_adjustment_factor
[1831]	380
[1320]	381	USE grid_variables, &
	382	ONLY: ddx2, ddy2
	383
	384	USE indices, &
[2232]	385	ONLY: nxlg, nxrg, nyng, nysg, nzb, nzb_max, nzt, wall_flags_0
[1320]	386
	387	USE kinds
[1849]	388
	389	USE microphysics_mod, &
	390	ONLY: collision_turbulence
	391
[1320]	392	USE particle_attributes, &
	393	ONLY: use_sgs_for_particles, wang_kernel
[1]	394
[2232]	395	USE surface_mod, &
	396	ONLY : bc_h
	397
[1]	398	IMPLICIT NONE
	399
[2232]	400	INTEGER(iwp) :: i !< running index x direction
	401	INTEGER(iwp) :: j !< running index y direction
	402	INTEGER(iwp) :: k !< running index z direction
	403	INTEGER(iwp) :: m !< running index surface elements
	404	INTEGER(iwp) :: surf_e !< End index of surface elements at (j,i)-gridpoint
	405	INTEGER(iwp) :: surf_s !< Start index of surface elements at (j,i)-gridpoint
	406
[1682]	407	REAL(wp) :: dvar_dz !<
[2232]	408	REAL(wp) :: flag !< flag to mask topography
[1682]	409	REAL(wp) :: l_stable !<
	410	REAL(wp) :: var_reference !<
[1]	411
[1010]	412	#if defined( __nopointer )
[1682]	413	REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: var !<
[1010]	414	#else
[1682]	415	REAL(wp), DIMENSION(:,:,:), POINTER :: var !<
[1010]	416	#endif
[1682]	417	REAL(wp), DIMENSION(nzb+1:nzt) :: dissipation !<
	418	REAL(wp), DIMENSION(nzb+1:nzt) :: l !<
	419	REAL(wp), DIMENSION(nzb+1:nzt) :: ll !<
[1]	420
	421	!
	422	!-- Calculate the mixing length (for dissipation)
[2232]	423	DO k = nzb+1, nzt
	424	!
	425	!-- Predetermine flag to mask topography
	426	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
	427
[97]	428	dvar_dz = atmos_ocean_sign * &
	429	( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
[1340]	430	IF ( dvar_dz > 0.0_wp ) THEN
[1179]	431	IF ( use_single_reference_value ) THEN
[1340]	432	l_stable = 0.76_wp * SQRT( e(k,j,i) ) / &
	433	SQRT( g / var_reference * dvar_dz ) + 1E-5_wp
[57]	434	ELSE
[1340]	435	l_stable = 0.76_wp * SQRT( e(k,j,i) ) / &
	436	SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5_wp
[57]	437	ENDIF
[1]	438	ELSE
	439	l_stable = l_grid(k)
	440	ENDIF
	441	!
	442	!-- Adjustment of the mixing length
	443	IF ( wall_adjustment ) THEN
[2232]	444	l(k) = MIN( wall_adjustment_factor * l_wall(k,j,i), &
	445	l_grid(k), l_stable )
	446	ll(k) = MIN( wall_adjustment_factor * l_wall(k,j,i), l_grid(k) )
[1]	447	ELSE
	448	l(k) = MIN( l_grid(k), l_stable )
	449	ll(k) = l_grid(k)
	450	ENDIF
	451	!
	452	!-- Calculate the tendency term
[2232]	453	dissipation(k) = ( 0.19_wp + 0.74_wp * l(k) / ll(k) ) * e(k,j,i) * &
[1340]	454	SQRT( e(k,j,i) ) / l(k)
[1]	455
[2232]	456	tend(k,j,i) = tend(k,j,i) &
	457	+ ( &
	458	( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) ) &
	459	- ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) ) &
	460	) * ddx2 * flag &
	461	+ ( &
	462	( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) ) &
	463	- ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) ) &
	464	) * ddy2 * flag &
	465	+ ( &
	466	( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
	467	* rho_air_zw(k) &
	468	- ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k) &
	469	* rho_air_zw(k-1) &
	470	) * ddzw(k) * drho_air(k) * flag &
	471	- dissipation(k) * flag
[1]	472
	473	ENDDO
	474
	475	!
	476	!-- Store dissipation if needed for calculating the sgs particle velocities
[1849]	477	IF ( use_sgs_for_particles .OR. wang_kernel .OR. &
	478	collision_turbulence ) THEN
[2232]	479	DO k = nzb+1, nzt
	480	diss(k,j,i) = dissipation(k) * &
	481	MERGE( 1.0_wp, 0.0_wp, &
	482	BTEST( wall_flags_0(k,j,i), 0 ) )
[1]	483	ENDDO
	484	!
[2232]	485	!-- Neumann boundary condition for dissipation diss(nzb,:,:) = diss(nzb+1,:,:)
	486	!-- For each surface type determine start and end index (in case of elevated
	487	!-- toopography several up/downward facing surfaces may exist.
	488	surf_s = bc_h(0)%start_index(j,i)
	489	surf_e = bc_h(0)%end_index(j,i)
	490	DO m = surf_s, surf_e
	491	k = bc_h(0)%k(m)
	492	diss(k-1,j,i) = diss(k,j,i)
	493	ENDDO
	494	!
	495	!-- Downward facing surfaces
	496	surf_s = bc_h(1)%start_index(j,i)
	497	surf_e = bc_h(1)%end_index(j,i)
	498	DO m = surf_s, surf_e
	499	k = bc_h(1)%k(m)
	500	diss(k+1,j,i) = diss(k,j,i)
	501	ENDDO
[1]	502	ENDIF
	503
	504	END SUBROUTINE diffusion_e_ij
	505
	506	END MODULE diffusion_e_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |