Home

Context Navigation

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

Last change on this file since 19 was 19, checked in by raasch, 18 years ago
preliminary version of modified boundary conditions at top
Property svn:keywords set to `Id`
File size: 9.1 KB

Rev	Line
[1]	1	MODULE diffusion_e_mod
	2
	3	!------------------------------------------------------------------------------!
	4	! Actual revisions:
	5	! -----------------
[19]	6	! Calculation extended for gridpoint nzt
[1]	7	!
	8	! Former revisions:
	9	! -----------------
[3]	10	! $Id: diffusion_e.f90 19 2007-02-23 04:53:48Z raasch $
	11	! RCS Log replace by Id keyword, revision history cleaned up
	12	!
[1]	13	! Revision 1.18 2006/08/04 14:29:43 raasch
	14	! dissipation is stored in extra array diss if needed later on for calculating
	15	! the sgs particle velocities
	16	!
	17	! Revision 1.1 1997/09/19 07:40:24 raasch
	18	! Initial revision
	19	!
	20	!
	21	! Description:
	22	! ------------
	23	! Diffusion- and dissipation terms for the TKE
	24	!------------------------------------------------------------------------------!
	25
	26	PRIVATE
	27	PUBLIC diffusion_e
	28
	29
	30	INTERFACE diffusion_e
	31	MODULE PROCEDURE diffusion_e
	32	MODULE PROCEDURE diffusion_e_ij
	33	END INTERFACE diffusion_e
	34
	35	CONTAINS
	36
	37
	38	!------------------------------------------------------------------------------!
	39	! Call for all grid points
	40	!------------------------------------------------------------------------------!
	41	SUBROUTINE diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, theta, &
	42	rif, tend, zu )
	43
	44	USE control_parameters
	45	USE grid_variables
	46	USE indices
	47	USE particle_attributes
	48
	49	IMPLICIT NONE
	50
	51	INTEGER :: i, j, k
	52	REAL :: dpt_dz, l_stable, phi_m
	53	REAL :: ddzu(1:nzt+1), dd2zu(1:nzt), ddzw(1:nzt), &
	54	l_grid(1:nzt), zu(0:nzt+1)
	55	REAL, DIMENSION(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) :: diss, tend
	56	REAL, DIMENSION(:,:), POINTER :: rif
	57	REAL, DIMENSION(:,:,:), POINTER :: e, km, theta
[19]	58	REAL, DIMENSION(nzb+1:nzt,nys:nyn) :: dissipation, l, ll
[1]	59
	60
	61	DO i = nxl, nxr
	62	DO j = nys, nyn
	63	!
	64	!-- First, calculate phi-function for eventually adjusting the &
	65	!-- mixing length to the prandtl mixing length
	66	IF ( adjust_mixing_length .AND. prandtl_layer ) THEN
	67	IF ( rif(j,i) >= 0.0 ) THEN
	68	phi_m = 1.0 + 5.0 * rif(j,i)
	69	ELSE
	70	phi_m = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rif(j,i) ) )
	71	ENDIF
	72	ENDIF
	73
[19]	74	DO k = nzb_s_inner(j,i)+1, nzt
[1]	75	!
	76	!-- Calculate the mixing length (for dissipation)
	77	dpt_dz = ( theta(k+1,j,i) - theta(k-1,j,i) ) * dd2zu(k)
	78	IF ( dpt_dz > 0.0 ) THEN
	79	l_stable = 0.76 * SQRT( e(k,j,i) ) / &
	80	SQRT( g / theta(k,j,i) * dpt_dz ) + 1E-5
	81	ELSE
	82	l_stable = l_grid(k)
	83	ENDIF
	84	!
	85	!-- Adjustment of the mixing length
	86	IF ( wall_adjustment ) THEN
	87	l(k,j) = MIN( wall_adjustment_factor * zu(k), l_grid(k), &
	88	l_stable )
	89	ll(k,j) = MIN( wall_adjustment_factor * zu(k), l_grid(k) )
	90	ELSE
	91	l(k,j) = MIN( l_grid(k), l_stable )
	92	ll(k,j) = l_grid(k)
	93	ENDIF
	94	IF ( adjust_mixing_length .AND. prandtl_layer ) THEN
	95	l(k,j) = MIN( l(k,j), kappa * zu(k) / phi_m )
	96	ll(k,j) = MIN( ll(k,j), kappa * zu(k) / phi_m )
	97	ENDIF
	98
	99	ENDDO
	100	ENDDO
	101	!
	102	!-- Calculate the tendency terms
	103	DO j = nys, nyn
[19]	104	DO k = nzb_s_inner(j,i)+1, nzt
[1]	105
	106	dissipation(k,j) = ( 0.19 + 0.74 * l(k,j) / ll(k,j) ) * &
	107	e(k,j,i) * SQRT( e(k,j,i) ) / l(k,j)
	108
	109	tend(k,j,i) = tend(k,j,i) &
	110	+ ( &
	111	( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) ) &
	112	- ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) ) &
	113	) * ddx2 &
	114	+ ( &
	115	( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) ) &
	116	- ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) ) &
	117	) * ddy2 &
	118	+ ( &
	119	( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
	120	- ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k) &
	121	) * ddzw(k) &
	122	- dissipation(k,j)
	123
	124	ENDDO
	125	ENDDO
	126
	127	!
	128	!-- Store dissipation if needed for calculating the sgs particle
	129	!-- velocities
	130	IF ( use_sgs_for_particles ) THEN
	131	DO j = nys, nyn
[19]	132	DO k = nzb_s_inner(j,i)+1, nzt
[1]	133	diss(k,j,i) = dissipation(k,j)
	134	ENDDO
	135	ENDDO
	136	ENDIF
	137
	138	ENDDO
	139
	140	!
	141	!-- Boundary condition for dissipation
	142	IF ( use_sgs_for_particles ) THEN
	143	DO i = nxl, nxr
	144	DO j = nys, nyn
	145	diss(nzb_s_inner(j,i),j,i) = diss(nzb_s_inner(j,i)+1,j,i)
	146	ENDDO
	147	ENDDO
	148	ENDIF
	149
	150	END SUBROUTINE diffusion_e
	151
	152
	153	!------------------------------------------------------------------------------!
	154	! Call for grid point i,j
	155	!------------------------------------------------------------------------------!
	156	SUBROUTINE diffusion_e_ij( i, j, ddzu, dd2zu, ddzw, diss, e, km, l_grid, &
	157	theta, rif, tend, zu )
	158
	159	USE control_parameters
	160	USE grid_variables
	161	USE indices
	162	USE particle_attributes
	163
	164	IMPLICIT NONE
	165
	166	INTEGER :: i, j, k
	167	REAL :: dpt_dz, l_stable, phi_m
	168	REAL :: ddzu(1:nzt+1), dd2zu(1:nzt), ddzw(1:nzt), &
	169	l_grid(1:nzt), zu(0:nzt+1)
	170	REAL, DIMENSION(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) :: diss, tend
	171	REAL, DIMENSION(:,:), POINTER :: rif
	172	REAL, DIMENSION(:,:,:), POINTER :: e, km, theta
[19]	173	REAL, DIMENSION(nzb+1:nzt) :: dissipation, l, ll
[1]	174
	175
	176	!
	177	!-- First, calculate phi-function for eventually adjusting the mixing length
	178	!-- to the prandtl mixing length
	179	IF ( adjust_mixing_length .AND. prandtl_layer ) THEN
	180	IF ( rif(j,i) >= 0.0 ) THEN
	181	phi_m = 1.0 + 5.0 * rif(j,i)
	182	ELSE
	183	phi_m = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rif(j,i) ) )
	184	ENDIF
	185	ENDIF
	186
	187	!
	188	!-- Calculate the mixing length (for dissipation)
[19]	189	DO k = nzb_s_inner(j,i)+1, nzt
[1]	190	dpt_dz = ( theta(k+1,j,i) - theta(k-1,j,i) ) * dd2zu(k)
	191	IF ( dpt_dz > 0.0 ) THEN
	192	l_stable = 0.76 * SQRT( e(k,j,i) ) / &
	193	SQRT( g / theta(k,j,i) * dpt_dz ) + 1E-5
	194	ELSE
	195	l_stable = l_grid(k)
	196	ENDIF
	197	!
	198	!-- Adjustment of the mixing length
	199	IF ( wall_adjustment ) THEN
	200	l(k) = MIN( wall_adjustment_factor * zu(k), l_grid(k), l_stable )
	201	ll(k) = MIN( wall_adjustment_factor * zu(k), l_grid(k) )
	202	ELSE
	203	l(k) = MIN( l_grid(k), l_stable )
	204	ll(k) = l_grid(k)
	205	ENDIF
	206	IF ( adjust_mixing_length .AND. prandtl_layer ) THEN
	207	l(k) = MIN( l(k), kappa * zu(k) / phi_m )
	208	ll(k) = MIN( ll(k), kappa * zu(k) / phi_m )
	209	ENDIF
	210
	211	!
	212	!-- Calculate the tendency term
	213	dissipation(k) = ( 0.19 + 0.74 * l(k) / ll(k) ) * e(k,j,i) * &
	214	SQRT( e(k,j,i) ) / l(k)
	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)
	230
	231	ENDDO
	232
	233	!
	234	!-- Store dissipation if needed for calculating the sgs particle velocities
	235	IF ( use_sgs_for_particles ) THEN
[19]	236	DO k = nzb_s_inner(j,i)+1, nzt
[1]	237	diss(k,j,i) = dissipation(k)
	238	ENDDO
	239	!
	240	!-- Boundary condition for dissipation
	241	diss(nzb_s_inner(j,i),j,i) = diss(nzb_s_inner(j,i)+1,j,i)
	242	ENDIF
	243
	244	END SUBROUTINE diffusion_e_ij
	245
	246	END MODULE diffusion_e_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |