Home

Context Navigation

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

Last change on this file since 1032 was 1017, checked in by raasch, 12 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 20.3 KB

Rev	Line
[1]	1	MODULE diffusion_e_mod
	2
	3	!------------------------------------------------------------------------------!
[484]	4	! Current revisions:
[1]	5	! -----------------
[98]	6	!
[1017]	7	!
[98]	8	! Former revisions:
	9	! -----------------
	10	! $Id: diffusion_e.f90 1017 2012-09-27 11:28:50Z letzel $
	11	!
[1017]	12	! 1015 2012-09-27 09:23:24Z raasch
	13	! accelerator version (*_acc) added,
	14	! adjustment of mixing length to the Prandtl mixing length at first grid point
	15	! above ground removed
	16	!
[1011]	17	! 1010 2012-09-20 07:59:54Z raasch
	18	! cpp switch __nopointer added for pointer free version
	19	!
[1002]	20	! 1001 2012-09-13 14:08:46Z raasch
	21	! most arrays comunicated by module instead of parameter list
	22	!
[826]	23	! 825 2012-02-19 03:03:44Z raasch
	24	! wang_collision_kernel renamed wang_kernel
	25	!
[791]	26	! 790 2011-11-29 03:11:20Z raasch
	27	! diss is also calculated in case that the Wang kernel is used
	28	!
[668]	29	! 667 2010-12-23 12:06:00Z suehring/gryschka
	30	! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng
	31	!
[98]	32	! 97 2007-06-21 08:23:15Z raasch
[94]	33	! Adjustment of mixing length calculation for the ocean version. zw added to
	34	! argument list.
	35	! This is also a bugfix, because the height above the topography is now
	36	! used instead of the height above level k=0.
[97]	37	! theta renamed var, dpt_dz renamed dvar_dz, +new argument var_reference
	38	! use_pt_reference renamed use_reference
[1]	39	!
[77]	40	! 65 2007-03-13 12:11:43Z raasch
	41	! Reference temperature pt_reference can be used in buoyancy term
	42	!
[39]	43	! 20 2007-02-26 00:12:32Z raasch
	44	! Bugfix: ddzw dimensioned 1:nzt"+1"
	45	! Calculation extended for gridpoint nzt
	46	!
[3]	47	! RCS Log replace by Id keyword, revision history cleaned up
	48	!
[1]	49	! Revision 1.18 2006/08/04 14:29:43 raasch
	50	! dissipation is stored in extra array diss if needed later on for calculating
	51	! the sgs particle velocities
	52	!
	53	! Revision 1.1 1997/09/19 07:40:24 raasch
	54	! Initial revision
	55	!
	56	!
	57	! Description:
	58	! ------------
	59	! Diffusion- and dissipation terms for the TKE
	60	!------------------------------------------------------------------------------!
	61
	62	PRIVATE
[1015]	63	PUBLIC diffusion_e, diffusion_e_acc
[1]	64
	65
	66	INTERFACE diffusion_e
	67	MODULE PROCEDURE diffusion_e
	68	MODULE PROCEDURE diffusion_e_ij
	69	END INTERFACE diffusion_e
	70
[1015]	71	INTERFACE diffusion_e_acc
	72	MODULE PROCEDURE diffusion_e_acc
	73	END INTERFACE diffusion_e_acc
	74
[1]	75	CONTAINS
	76
	77
	78	!------------------------------------------------------------------------------!
	79	! Call for all grid points
	80	!------------------------------------------------------------------------------!
[1001]	81	SUBROUTINE diffusion_e( var, var_reference )
[1]	82
[1001]	83	USE arrays_3d
[1]	84	USE control_parameters
	85	USE grid_variables
	86	USE indices
	87	USE particle_attributes
	88
	89	IMPLICIT NONE
	90
	91	INTEGER :: i, j, k
[1015]	92	REAL :: dvar_dz, l_stable, var_reference
[1001]	93
[1010]	94	#if defined( __nopointer )
	95	REAL, DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: var
	96	#else
[1001]	97	REAL, DIMENSION(:,:,:), POINTER :: var
[1010]	98	#endif
[19]	99	REAL, DIMENSION(nzb+1:nzt,nys:nyn) :: dissipation, l, ll
[1]	100
	101
	102	!
[65]	103	!-- This if clause must be outside the k-loop because otherwise
	104	!-- runtime errors occur with -C hopt on NEC
[97]	105	IF ( use_reference ) THEN
[65]	106
	107	DO i = nxl, nxr
	108	DO j = nys, nyn
	109	DO k = nzb_s_inner(j,i)+1, nzt
[1]	110	!
[65]	111	!-- Calculate the mixing length (for dissipation)
[97]	112	dvar_dz = atmos_ocean_sign * &
	113	( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
	114	IF ( dvar_dz > 0.0 ) THEN
[57]	115	l_stable = 0.76 * SQRT( e(k,j,i) ) / &
[97]	116	SQRT( g / var_reference * dvar_dz ) + 1E-5
[57]	117	ELSE
[65]	118	l_stable = l_grid(k)
[57]	119	ENDIF
[1]	120	!
[65]	121	!-- Adjustment of the mixing length
	122	IF ( wall_adjustment ) THEN
[94]	123	l(k,j) = MIN( wall_adjustment_factor * &
	124	( zu(k) - zw(nzb_s_inner(j,i)) ), &
	125	l_grid(k), l_stable )
	126	ll(k,j) = MIN( wall_adjustment_factor * &
	127	( zu(k) - zw(nzb_s_inner(j,i)) ), &
	128	l_grid(k) )
[65]	129	ELSE
	130	l(k,j) = MIN( l_grid(k), l_stable )
	131	ll(k,j) = l_grid(k)
	132	ENDIF
[1]	133
[65]	134	ENDDO
[1]	135	ENDDO
[65]	136
[1]	137	!
[65]	138	!-- Calculate the tendency terms
	139	DO j = nys, nyn
	140	DO k = nzb_s_inner(j,i)+1, nzt
[1]	141
[65]	142	dissipation(k,j) = ( 0.19 + 0.74 * l(k,j) / ll(k,j) ) * &
	143	e(k,j,i) * SQRT( e(k,j,i) ) / l(k,j)
[1]	144
[65]	145	tend(k,j,i) = tend(k,j,i) &
[1]	146	+ ( &
	147	( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) ) &
	148	- ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) ) &
	149	) * ddx2 &
	150	+ ( &
	151	( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) ) &
	152	- ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) ) &
	153	) * ddy2 &
	154	+ ( &
	155	( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
	156	- ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k) &
	157	) * ddzw(k) &
	158	- dissipation(k,j)
	159
[65]	160	ENDDO
[1]	161	ENDDO
[65]	162
	163	!
	164	!-- Store dissipation if needed for calculating the sgs particle
	165	!-- velocities
[825]	166	IF ( use_sgs_for_particles .OR. wang_kernel ) THEN
[65]	167	DO j = nys, nyn
	168	DO k = nzb_s_inner(j,i)+1, nzt
	169	diss(k,j,i) = dissipation(k,j)
	170	ENDDO
	171	ENDDO
	172	ENDIF
	173
[1]	174	ENDDO
	175
[65]	176	ELSE
	177
	178	DO i = nxl, nxr
	179	DO j = nys, nyn
	180	DO k = nzb_s_inner(j,i)+1, nzt
	181	!
	182	!-- Calculate the mixing length (for dissipation)
[97]	183	dvar_dz = atmos_ocean_sign * &
	184	( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
	185	IF ( dvar_dz > 0.0 ) THEN
[65]	186	l_stable = 0.76 * SQRT( e(k,j,i) ) / &
[97]	187	SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5
[65]	188	ELSE
	189	l_stable = l_grid(k)
	190	ENDIF
	191	!
	192	!-- Adjustment of the mixing length
	193	IF ( wall_adjustment ) THEN
[94]	194	l(k,j) = MIN( wall_adjustment_factor * &
	195	( zu(k) - zw(nzb_s_inner(j,i)) ), &
	196	l_grid(k), l_stable )
	197	ll(k,j) = MIN( wall_adjustment_factor * &
	198	( zu(k) - zw(nzb_s_inner(j,i)) ), &
	199	l_grid(k) )
[65]	200	ELSE
	201	l(k,j) = MIN( l_grid(k), l_stable )
	202	ll(k,j) = l_grid(k)
	203	ENDIF
	204
	205	ENDDO
	206	ENDDO
	207
	208	!
	209	!-- Calculate the tendency terms
[1]	210	DO j = nys, nyn
[19]	211	DO k = nzb_s_inner(j,i)+1, nzt
[65]	212
	213	dissipation(k,j) = ( 0.19 + 0.74 * l(k,j) / ll(k,j) ) * &
	214	e(k,j,i) * SQRT( e(k,j,i) ) / l(k,j)
	215
	216	tend(k,j,i) = tend(k,j,i) &
	217	+ ( &
	218	( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) ) &
	219	- ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) ) &
	220	) * ddx2 &
	221	+ ( &
	222	( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) ) &
	223	- ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) ) &
	224	) * ddy2 &
	225	+ ( &
	226	( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
	227	- ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k) &
	228	) * ddzw(k) &
	229	- dissipation(k,j)
	230
[1]	231	ENDDO
	232	ENDDO
	233
[65]	234	!
	235	!-- Store dissipation if needed for calculating the sgs particle
	236	!-- velocities
[825]	237	IF ( use_sgs_for_particles .OR. wang_kernel ) THEN
[65]	238	DO j = nys, nyn
	239	DO k = nzb_s_inner(j,i)+1, nzt
	240	diss(k,j,i) = dissipation(k,j)
	241	ENDDO
	242	ENDDO
	243	ENDIF
[1]	244
[65]	245	ENDDO
	246
	247	ENDIF
	248
[1]	249	!
	250	!-- Boundary condition for dissipation
[825]	251	IF ( use_sgs_for_particles .OR. wang_kernel ) THEN
[1]	252	DO i = nxl, nxr
	253	DO j = nys, nyn
	254	diss(nzb_s_inner(j,i),j,i) = diss(nzb_s_inner(j,i)+1,j,i)
	255	ENDDO
	256	ENDDO
	257	ENDIF
	258
	259	END SUBROUTINE diffusion_e
	260
	261
	262	!------------------------------------------------------------------------------!
[1015]	263	! Call for all grid points - accelerator version
	264	!------------------------------------------------------------------------------!
	265	SUBROUTINE diffusion_e_acc( var, var_reference )
	266
	267	USE arrays_3d
	268	USE control_parameters
	269	USE grid_variables
	270	USE indices
	271	USE particle_attributes
	272
	273	IMPLICIT NONE
	274
	275	INTEGER :: i, j, k
	276	REAL :: dissipation, dvar_dz, l, ll, l_stable, var_reference
	277
	278	#if defined( __nopointer )
	279	REAL, DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: var
	280	#else
	281	REAL, DIMENSION(:,:,:), POINTER :: var
	282	#endif
	283
	284
	285	!
	286	!-- This if clause must be outside the k-loop because otherwise
	287	!-- runtime errors occur with -C hopt on NEC
	288	IF ( use_reference ) THEN
	289	STOP '+++ use_reference in diffusion_e not implemented'
	290	! DO i = nxl, nxr
	291	! DO j = nys, nyn
	292	! DO k = nzb_s_inner(j,i)+1, nzt
	293	!
	294	!-- Calculate the mixing length (for dissipation)
	295	! dvar_dz = atmos_ocean_sign * &
	296	! ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
	297	! IF ( dvar_dz > 0.0 ) THEN
	298	! l_stable = 0.76 * SQRT( e(k,j,i) ) / &
	299	! SQRT( g / var_reference * dvar_dz ) + 1E-5
	300	! ELSE
	301	! l_stable = l_grid(k)
	302	! ENDIF
	303	!
	304	!-- Adjustment of the mixing length
	305	! IF ( wall_adjustment ) THEN
	306	! l(k,j) = MIN( wall_adjustment_factor * &
	307	! ( zu(k) - zw(nzb_s_inner(j,i)) ), &
	308	! l_grid(k), l_stable )
	309	! ll(k,j) = MIN( wall_adjustment_factor * &
	310	! ( zu(k) - zw(nzb_s_inner(j,i)) ), &
	311	! l_grid(k) )
	312	! ELSE
	313	! l(k,j) = MIN( l_grid(k), l_stable )
	314	! ll(k,j) = l_grid(k)
	315	! ENDIF
	316	!
	317	! ENDDO
	318	! ENDDO
	319	!
	320	!
	321	!-- Calculate the tendency terms
	322	! DO j = nys, nyn
	323	! DO k = nzb_s_inner(j,i)+1, nzt
	324	!
	325	! dissipation(k,j) = ( 0.19 + 0.74 * l(k,j) / ll(k,j) ) * &
	326	! e(k,j,i) * SQRT( e(k,j,i) ) / l(k,j)
	327	!
	328	! tend(k,j,i) = tend(k,j,i) &
	329	! + ( &
	330	! ( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) ) &
	331	! - ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) ) &
	332	! ) * ddx2 &
	333	! + ( &
	334	! ( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) ) &
	335	! - ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) ) &
	336	! ) * ddy2 &
	337	! + ( &
	338	! ( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
	339	! - ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k) &
	340	! ) * ddzw(k) &
	341	! - dissipation(k,j)
	342	!
	343	! ENDDO
	344	! ENDDO
	345	!
	346	!
	347	!-- Store dissipation if needed for calculating the sgs particle
	348	!-- velocities
	349	! IF ( use_sgs_for_particles .OR. wang_kernel ) THEN
	350	! DO j = nys, nyn
	351	! DO k = nzb_s_inner(j,i)+1, nzt
	352	! diss(k,j,i) = dissipation(k,j)
	353	! ENDDO
	354	! ENDDO
	355	! ENDIF
	356	!
	357	! ENDDO
	358	!
	359	ELSE
	360
	361	!$acc kernels present( ddzu, ddzw, dd2zu, diss, e, km, l_grid ) &
	362	!$acc present( nzb_s_inner, rif, tend, var, zu, zw )
	363	!$acc loop
	364	DO i = nxl, nxr
	365	DO j = nys, nyn
	366	!$acc loop vector( 32 )
	367	DO k = 1, nzt
	368
	369	IF ( k > nzb_s_inner(j,i) ) THEN
	370	!
	371	!-- Calculate the mixing length (for dissipation)
	372	dvar_dz = atmos_ocean_sign * &
	373	( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
	374	IF ( dvar_dz > 0.0 ) THEN
	375	l_stable = 0.76 * SQRT( e(k,j,i) ) / &
	376	SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5
	377	ELSE
	378	l_stable = l_grid(k)
	379	ENDIF
	380	!
	381	!-- Adjustment of the mixing length
	382	IF ( wall_adjustment ) THEN
	383	l = MIN( wall_adjustment_factor * &
	384	( zu(k) - zw(nzb_s_inner(j,i)) ), &
	385	l_grid(k), l_stable )
	386	ll = MIN( wall_adjustment_factor * &
	387	( zu(k) - zw(nzb_s_inner(j,i)) ), &
	388	l_grid(k) )
	389	ELSE
	390	l = MIN( l_grid(k), l_stable )
	391	ll = l_grid(k)
	392	ENDIF
	393	!
	394	!-- Calculate the tendency terms
	395	dissipation = ( 0.19 + 0.74 * l / ll ) * &
	396	e(k,j,i) * SQRT( e(k,j,i) ) / l
	397
	398	tend(k,j,i) = tend(k,j,i) &
	399	+ ( &
	400	( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) ) &
	401	- ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) ) &
	402	) * ddx2 &
	403	+ ( &
	404	( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) ) &
	405	- ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) ) &
	406	) * ddy2 &
	407	+ ( &
	408	( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
	409	- ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k) &
	410	) * ddzw(k) &
	411	- dissipation
	412
	413	!
	414	!-- Store dissipation if needed for calculating the sgs
	415	!-- particle velocities
	416	IF ( use_sgs_for_particles .OR. wang_kernel ) THEN
	417	diss(k,j,i) = dissipation
	418	ENDIF
	419
	420	ENDIF
	421
	422	ENDDO
	423	ENDDO
	424	ENDDO
	425	!$acc end kernels
	426
	427	ENDIF
	428
	429	!
	430	!-- Boundary condition for dissipation
	431	IF ( use_sgs_for_particles .OR. wang_kernel ) THEN
	432	!$acc kernels present( diss, nzb_s_inner )
	433	!$acc loop
	434	DO i = nxl, nxr
	435	!$acc loop vector( 32 )
	436	DO j = nys, nyn
	437	diss(nzb_s_inner(j,i),j,i) = diss(nzb_s_inner(j,i)+1,j,i)
	438	ENDDO
	439	ENDDO
	440	!$acc end kernels
	441	ENDIF
	442
	443	END SUBROUTINE diffusion_e_acc
	444
	445
	446	!------------------------------------------------------------------------------!
[1]	447	! Call for grid point i,j
	448	!------------------------------------------------------------------------------!
[1001]	449	SUBROUTINE diffusion_e_ij( i, j, var, var_reference )
[1]	450
[1001]	451	USE arrays_3d
[1]	452	USE control_parameters
	453	USE grid_variables
	454	USE indices
	455	USE particle_attributes
	456
	457	IMPLICIT NONE
	458
[1001]	459	INTEGER :: i, j, k
[1015]	460	REAL :: dvar_dz, l_stable, var_reference
[1]	461
[1010]	462	#if defined( __nopointer )
	463	REAL, DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: var
	464	#else
[1001]	465	REAL, DIMENSION(:,:,:), POINTER :: var
[1010]	466	#endif
[1001]	467	REAL, DIMENSION(nzb+1:nzt) :: dissipation, l, ll
[1]	468
[1001]	469
[1]	470	!
	471	!-- Calculate the mixing length (for dissipation)
[19]	472	DO k = nzb_s_inner(j,i)+1, nzt
[97]	473	dvar_dz = atmos_ocean_sign * &
	474	( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
	475	IF ( dvar_dz > 0.0 ) THEN
	476	IF ( use_reference ) THEN
[57]	477	l_stable = 0.76 * SQRT( e(k,j,i) ) / &
[97]	478	SQRT( g / var_reference * dvar_dz ) + 1E-5
[57]	479	ELSE
	480	l_stable = 0.76 * SQRT( e(k,j,i) ) / &
[97]	481	SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5
[57]	482	ENDIF
[1]	483	ELSE
	484	l_stable = l_grid(k)
	485	ENDIF
	486	!
	487	!-- Adjustment of the mixing length
	488	IF ( wall_adjustment ) THEN
[94]	489	l(k) = MIN( wall_adjustment_factor * &
	490	( zu(k) - zw(nzb_s_inner(j,i)) ), l_grid(k), &
	491	l_stable )
	492	ll(k) = MIN( wall_adjustment_factor * &
	493	( zu(k) - zw(nzb_s_inner(j,i)) ), l_grid(k) )
[1]	494	ELSE
	495	l(k) = MIN( l_grid(k), l_stable )
	496	ll(k) = l_grid(k)
	497	ENDIF
	498	!
	499	!-- Calculate the tendency term
	500	dissipation(k) = ( 0.19 + 0.74 * l(k) / ll(k) ) * e(k,j,i) * &
	501	SQRT( e(k,j,i) ) / l(k)
	502
	503	tend(k,j,i) = tend(k,j,i) &
	504	+ ( &
	505	( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) ) &
	506	- ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) ) &
	507	) * ddx2 &
	508	+ ( &
	509	( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) ) &
	510	- ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) ) &
	511	) * ddy2 &
	512	+ ( &
	513	( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
	514	- ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k) &
	515	) * ddzw(k) &
	516	- dissipation(k)
	517
	518	ENDDO
	519
	520	!
	521	!-- Store dissipation if needed for calculating the sgs particle velocities
[825]	522	IF ( use_sgs_for_particles .OR. wang_kernel ) THEN
[19]	523	DO k = nzb_s_inner(j,i)+1, nzt
[1]	524	diss(k,j,i) = dissipation(k)
	525	ENDDO
	526	!
	527	!-- Boundary condition for dissipation
	528	diss(nzb_s_inner(j,i),j,i) = diss(nzb_s_inner(j,i)+1,j,i)
	529	ENDIF
	530
	531	END SUBROUTINE diffusion_e_ij
	532
	533	END MODULE diffusion_e_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |