Home

Context Navigation

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

Last change on this file since 1004 was 1002, checked in by raasch, 12 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 14.8 KB

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