Home

Context Navigation

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

Last change on this file since 738 was 668, checked in by suehring, 14 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 15.1 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |