Home

Context Navigation

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

Last change on this file since 1016 was 1015, checked in by raasch, 12 years ago
Starting with changes required for GPU optimization. OpenACC statements for using NVIDIA GPUs added. Adjustment of mixing length to the Prandtl mixing length at first grid point above ground removed. mask array is set zero for ghost boundaries
Property svn:keywords set to `Id`
File size: 20.2 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |