Home

Context Navigation

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

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |