Home

Context Navigation

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

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

Rev	Line
[1]	1	MODULE diffusion_e_mod
	2
	3	!------------------------------------------------------------------------------!
[484]	4	! Current revisions:
[1]	5	! -----------------
[98]	6	!
[1011]	7	!
[98]	8	! Former revisions:
	9	! -----------------
	10	! $Id: diffusion_e.f90 1011 2012-09-20 08:15:29Z raasch $
	11	!
[1011]	12	! 1010 2012-09-20 07:59:54Z raasch
	13	! cpp switch __nopointer added for pointer free version
	14	!
[1002]	15	! 1001 2012-09-13 14:08:46Z raasch
	16	! most arrays comunicated by module instead of parameter list
	17	!
[826]	18	! 825 2012-02-19 03:03:44Z raasch
	19	! wang_collision_kernel renamed wang_kernel
	20	!
[791]	21	! 790 2011-11-29 03:11:20Z raasch
	22	! diss is also calculated in case that the Wang kernel is used
	23	!
[668]	24	! 667 2010-12-23 12:06:00Z suehring/gryschka
	25	! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng
	26	!
[98]	27	! 97 2007-06-21 08:23:15Z raasch
[94]	28	! Adjustment of mixing length calculation for the ocean version. zw added to
	29	! argument list.
	30	! This is also a bugfix, because the height above the topography is now
	31	! used instead of the height above level k=0.
[97]	32	! theta renamed var, dpt_dz renamed dvar_dz, +new argument var_reference
	33	! use_pt_reference renamed use_reference
[1]	34	!
[77]	35	! 65 2007-03-13 12:11:43Z raasch
	36	! Reference temperature pt_reference can be used in buoyancy term
	37	!
[39]	38	! 20 2007-02-26 00:12:32Z raasch
	39	! Bugfix: ddzw dimensioned 1:nzt"+1"
	40	! Calculation extended for gridpoint nzt
	41	!
[3]	42	! RCS Log replace by Id keyword, revision history cleaned up
	43	!
[1]	44	! Revision 1.18 2006/08/04 14:29:43 raasch
	45	! dissipation is stored in extra array diss if needed later on for calculating
	46	! the sgs particle velocities
	47	!
	48	! Revision 1.1 1997/09/19 07:40:24 raasch
	49	! Initial revision
	50	!
	51	!
	52	! Description:
	53	! ------------
	54	! Diffusion- and dissipation terms for the TKE
	55	!------------------------------------------------------------------------------!
	56
	57	PRIVATE
	58	PUBLIC diffusion_e
	59
	60
	61	INTERFACE diffusion_e
	62	MODULE PROCEDURE diffusion_e
	63	MODULE PROCEDURE diffusion_e_ij
	64	END INTERFACE diffusion_e
	65
	66	CONTAINS
	67
	68
	69	!------------------------------------------------------------------------------!
	70	! Call for all grid points
	71	!------------------------------------------------------------------------------!
[1001]	72	SUBROUTINE diffusion_e( var, var_reference )
[1]	73
[1001]	74	USE arrays_3d
[1]	75	USE control_parameters
	76	USE grid_variables
	77	USE indices
	78	USE particle_attributes
	79
	80	IMPLICIT NONE
	81
	82	INTEGER :: i, j, k
[1001]	83	REAL :: dvar_dz, l_stable, phi_m, var_reference
	84
[1010]	85	#if defined( __nopointer )
	86	REAL, DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: var
	87	#else
[1001]	88	REAL, DIMENSION(:,:,:), POINTER :: var
[1010]	89	#endif
[19]	90	REAL, DIMENSION(nzb+1:nzt,nys:nyn) :: dissipation, l, ll
[1]	91
	92
	93	!
[65]	94	!-- This if clause must be outside the k-loop because otherwise
	95	!-- runtime errors occur with -C hopt on NEC
[97]	96	IF ( use_reference ) THEN
[65]	97
	98	DO i = nxl, nxr
	99	DO j = nys, nyn
	100	!
	101	!-- First, calculate phi-function for eventually adjusting the &
	102	!-- mixing length to the prandtl mixing length
	103	IF ( adjust_mixing_length .AND. prandtl_layer ) THEN
	104	IF ( rif(j,i) >= 0.0 ) THEN
	105	phi_m = 1.0 + 5.0 * rif(j,i)
	106	ELSE
	107	phi_m = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rif(j,i) ) )
	108	ENDIF
[1]	109	ENDIF
	110
[65]	111	DO k = nzb_s_inner(j,i)+1, nzt
[1]	112	!
[65]	113	!-- Calculate the mixing length (for dissipation)
[97]	114	dvar_dz = atmos_ocean_sign * &
	115	( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
	116	IF ( dvar_dz > 0.0 ) THEN
[57]	117	l_stable = 0.76 * SQRT( e(k,j,i) ) / &
[97]	118	SQRT( g / var_reference * dvar_dz ) + 1E-5
[57]	119	ELSE
[65]	120	l_stable = l_grid(k)
[57]	121	ENDIF
[1]	122	!
[65]	123	!-- Adjustment of the mixing length
	124	IF ( wall_adjustment ) THEN
[94]	125	l(k,j) = MIN( wall_adjustment_factor * &
	126	( zu(k) - zw(nzb_s_inner(j,i)) ), &
	127	l_grid(k), l_stable )
	128	ll(k,j) = MIN( wall_adjustment_factor * &
	129	( zu(k) - zw(nzb_s_inner(j,i)) ), &
	130	l_grid(k) )
[65]	131	ELSE
	132	l(k,j) = MIN( l_grid(k), l_stable )
	133	ll(k,j) = l_grid(k)
	134	ENDIF
	135	IF ( adjust_mixing_length .AND. prandtl_layer ) THEN
[94]	136	l(k,j) = MIN( l(k,j), kappa * &
	137	( zu(k) - zw(nzb_s_inner(j,i)) ) &
	138	/ phi_m )
	139	ll(k,j) = MIN( ll(k,j), kappa * &
	140	( zu(k) - zw(nzb_s_inner(j,i)) ) &
	141	/ phi_m )
[65]	142	ENDIF
[1]	143
[65]	144	ENDDO
[1]	145	ENDDO
[65]	146
[1]	147	!
[65]	148	!-- Calculate the tendency terms
	149	DO j = nys, nyn
	150	DO k = nzb_s_inner(j,i)+1, nzt
[1]	151
[65]	152	dissipation(k,j) = ( 0.19 + 0.74 * l(k,j) / ll(k,j) ) * &
	153	e(k,j,i) * SQRT( e(k,j,i) ) / l(k,j)
[1]	154
[65]	155	tend(k,j,i) = tend(k,j,i) &
[1]	156	+ ( &
	157	( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) ) &
	158	- ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) ) &
	159	) * ddx2 &
	160	+ ( &
	161	( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) ) &
	162	- ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) ) &
	163	) * ddy2 &
	164	+ ( &
	165	( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
	166	- ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k) &
	167	) * ddzw(k) &
	168	- dissipation(k,j)
	169
[65]	170	ENDDO
[1]	171	ENDDO
[65]	172
	173	!
	174	!-- Store dissipation if needed for calculating the sgs particle
	175	!-- velocities
[825]	176	IF ( use_sgs_for_particles .OR. wang_kernel ) THEN
[65]	177	DO j = nys, nyn
	178	DO k = nzb_s_inner(j,i)+1, nzt
	179	diss(k,j,i) = dissipation(k,j)
	180	ENDDO
	181	ENDDO
	182	ENDIF
	183
[1]	184	ENDDO
	185
[65]	186	ELSE
	187
	188	DO i = nxl, nxr
	189	DO j = nys, nyn
[1]	190	!
[65]	191	!-- First, calculate phi-function for eventually adjusting the &
	192	!-- mixing length to the prandtl mixing length
	193	IF ( adjust_mixing_length .AND. prandtl_layer ) THEN
	194	IF ( rif(j,i) >= 0.0 ) THEN
	195	phi_m = 1.0 + 5.0 * rif(j,i)
	196	ELSE
	197	phi_m = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rif(j,i) ) )
	198	ENDIF
	199	ENDIF
	200
	201	DO k = nzb_s_inner(j,i)+1, nzt
	202	!
	203	!-- Calculate the mixing length (for dissipation)
[97]	204	dvar_dz = atmos_ocean_sign * &
	205	( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
	206	IF ( dvar_dz > 0.0 ) THEN
[65]	207	l_stable = 0.76 * SQRT( e(k,j,i) ) / &
[97]	208	SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5
[65]	209	ELSE
	210	l_stable = l_grid(k)
	211	ENDIF
	212	!
	213	!-- Adjustment of the mixing length
	214	IF ( wall_adjustment ) THEN
[94]	215	l(k,j) = MIN( wall_adjustment_factor * &
	216	( zu(k) - zw(nzb_s_inner(j,i)) ), &
	217	l_grid(k), l_stable )
	218	ll(k,j) = MIN( wall_adjustment_factor * &
	219	( zu(k) - zw(nzb_s_inner(j,i)) ), &
	220	l_grid(k) )
[65]	221	ELSE
	222	l(k,j) = MIN( l_grid(k), l_stable )
	223	ll(k,j) = l_grid(k)
	224	ENDIF
	225	IF ( adjust_mixing_length .AND. prandtl_layer ) THEN
[94]	226	l(k,j) = MIN( l(k,j), kappa * &
	227	( zu(k) - zw(nzb_s_inner(j,i)) ) &
	228	/ phi_m )
	229	ll(k,j) = MIN( ll(k,j), kappa * &
	230	( zu(k) - zw(nzb_s_inner(j,i)) ) &
	231	/ phi_m )
[65]	232	ENDIF
	233
	234	ENDDO
	235	ENDDO
	236
	237	!
	238	!-- Calculate the tendency terms
[1]	239	DO j = nys, nyn
[19]	240	DO k = nzb_s_inner(j,i)+1, nzt
[65]	241
	242	dissipation(k,j) = ( 0.19 + 0.74 * l(k,j) / ll(k,j) ) * &
	243	e(k,j,i) * SQRT( e(k,j,i) ) / l(k,j)
	244
	245	tend(k,j,i) = tend(k,j,i) &
	246	+ ( &
	247	( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) ) &
	248	- ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) ) &
	249	) * ddx2 &
	250	+ ( &
	251	( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) ) &
	252	- ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) ) &
	253	) * ddy2 &
	254	+ ( &
	255	( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
	256	- ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k) &
	257	) * ddzw(k) &
	258	- dissipation(k,j)
	259
[1]	260	ENDDO
	261	ENDDO
	262
[65]	263	!
	264	!-- Store dissipation if needed for calculating the sgs particle
	265	!-- velocities
[825]	266	IF ( use_sgs_for_particles .OR. wang_kernel ) THEN
[65]	267	DO j = nys, nyn
	268	DO k = nzb_s_inner(j,i)+1, nzt
	269	diss(k,j,i) = dissipation(k,j)
	270	ENDDO
	271	ENDDO
	272	ENDIF
[1]	273
[65]	274	ENDDO
	275
	276	ENDIF
	277
[1]	278	!
	279	!-- Boundary condition for dissipation
[825]	280	IF ( use_sgs_for_particles .OR. wang_kernel ) THEN
[1]	281	DO i = nxl, nxr
	282	DO j = nys, nyn
	283	diss(nzb_s_inner(j,i),j,i) = diss(nzb_s_inner(j,i)+1,j,i)
	284	ENDDO
	285	ENDDO
	286	ENDIF
	287
	288	END SUBROUTINE diffusion_e
	289
	290
	291	!------------------------------------------------------------------------------!
	292	! Call for grid point i,j
	293	!------------------------------------------------------------------------------!
[1001]	294	SUBROUTINE diffusion_e_ij( i, j, var, var_reference )
[1]	295
[1001]	296	USE arrays_3d
[1]	297	USE control_parameters
	298	USE grid_variables
	299	USE indices
	300	USE particle_attributes
	301
	302	IMPLICIT NONE
	303
[1001]	304	INTEGER :: i, j, k
	305	REAL :: dvar_dz, l_stable, phi_m, var_reference
[1]	306
[1010]	307	#if defined( __nopointer )
	308	REAL, DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: var
	309	#else
[1001]	310	REAL, DIMENSION(:,:,:), POINTER :: var
[1010]	311	#endif
[1001]	312	REAL, DIMENSION(nzb+1:nzt) :: dissipation, l, ll
[1]	313
[1001]	314
[1]	315	!
	316	!-- First, calculate phi-function for eventually adjusting the mixing length
	317	!-- to the prandtl mixing length
	318	IF ( adjust_mixing_length .AND. prandtl_layer ) THEN
	319	IF ( rif(j,i) >= 0.0 ) THEN
	320	phi_m = 1.0 + 5.0 * rif(j,i)
	321	ELSE
	322	phi_m = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rif(j,i) ) )
	323	ENDIF
	324	ENDIF
	325
	326	!
	327	!-- Calculate the mixing length (for dissipation)
[19]	328	DO k = nzb_s_inner(j,i)+1, nzt
[97]	329	dvar_dz = atmos_ocean_sign * &
	330	( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
	331	IF ( dvar_dz > 0.0 ) THEN
	332	IF ( use_reference ) THEN
[57]	333	l_stable = 0.76 * SQRT( e(k,j,i) ) / &
[97]	334	SQRT( g / var_reference * dvar_dz ) + 1E-5
[57]	335	ELSE
	336	l_stable = 0.76 * SQRT( e(k,j,i) ) / &
[97]	337	SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5
[57]	338	ENDIF
[1]	339	ELSE
	340	l_stable = l_grid(k)
	341	ENDIF
	342	!
	343	!-- Adjustment of the mixing length
	344	IF ( wall_adjustment ) THEN
[94]	345	l(k) = MIN( wall_adjustment_factor * &
	346	( zu(k) - zw(nzb_s_inner(j,i)) ), l_grid(k), &
	347	l_stable )
	348	ll(k) = MIN( wall_adjustment_factor * &
	349	( zu(k) - zw(nzb_s_inner(j,i)) ), l_grid(k) )
[1]	350	ELSE
	351	l(k) = MIN( l_grid(k), l_stable )
	352	ll(k) = l_grid(k)
	353	ENDIF
	354	IF ( adjust_mixing_length .AND. prandtl_layer ) THEN
[94]	355	l(k) = MIN( l(k), kappa * &
	356	( zu(k) - zw(nzb_s_inner(j,i)) ) / phi_m )
	357	ll(k) = MIN( ll(k), kappa * &
	358	( zu(k) - zw(nzb_s_inner(j,i)) ) / phi_m )
[1]	359	ENDIF
	360
	361	!
	362	!-- Calculate the tendency term
	363	dissipation(k) = ( 0.19 + 0.74 * l(k) / ll(k) ) * e(k,j,i) * &
	364	SQRT( e(k,j,i) ) / l(k)
	365
	366	tend(k,j,i) = tend(k,j,i) &
	367	+ ( &
	368	( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) ) &
	369	- ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) ) &
	370	) * ddx2 &
	371	+ ( &
	372	( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) ) &
	373	- ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) ) &
	374	) * ddy2 &
	375	+ ( &
	376	( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
	377	- ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k) &
	378	) * ddzw(k) &
	379	- dissipation(k)
	380
	381	ENDDO
	382
	383	!
	384	!-- Store dissipation if needed for calculating the sgs particle velocities
[825]	385	IF ( use_sgs_for_particles .OR. wang_kernel ) THEN
[19]	386	DO k = nzb_s_inner(j,i)+1, nzt
[1]	387	diss(k,j,i) = dissipation(k)
	388	ENDDO
	389	!
	390	!-- Boundary condition for dissipation
	391	diss(nzb_s_inner(j,i),j,i) = diss(nzb_s_inner(j,i)+1,j,i)
	392	ENDIF
	393
	394	END SUBROUTINE diffusion_e_ij
	395
	396	END MODULE diffusion_e_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |