Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |