Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |