Home

Context Navigation

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

Last change on this file since 513 was 484, checked in by raasch, 15 years ago
typo in file headers removed
Property svn:keywords set to `Id`
File size: 15.0 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |