Home

Context Navigation

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

Last change on this file since 1010 was 1010, checked in by raasch, 12 years ago
pointer free version can be generated with cpp switch nopointer
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	! -----------------
[1010]	6	! cpp switch __nopointer added for pointer free version
[98]	7	!
	8	! Former revisions:
	9	! -----------------
	10	! $Id: diffusion_e.f90 1010 2012-09-20 07:59:54Z 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
[1010]	82	#if defined( __nopointer )
	83	REAL, DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: var
	84	#else
[1001]	85	REAL, DIMENSION(:,:,:), POINTER :: var
[1010]	86	#endif
[19]	87	REAL, DIMENSION(nzb+1:nzt,nys:nyn) :: dissipation, l, ll
[1]	88
	89
	90	!
[65]	91	!-- This if clause must be outside the k-loop because otherwise
	92	!-- runtime errors occur with -C hopt on NEC
[97]	93	IF ( use_reference ) THEN
[65]	94
	95	DO i = nxl, nxr
	96	DO j = nys, nyn
	97	!
	98	!-- First, calculate phi-function for eventually adjusting the &
	99	!-- mixing length to the prandtl mixing length
	100	IF ( adjust_mixing_length .AND. prandtl_layer ) THEN
	101	IF ( rif(j,i) >= 0.0 ) THEN
	102	phi_m = 1.0 + 5.0 * rif(j,i)
	103	ELSE
	104	phi_m = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rif(j,i) ) )
	105	ENDIF
[1]	106	ENDIF
	107
[65]	108	DO k = nzb_s_inner(j,i)+1, nzt
[1]	109	!
[65]	110	!-- Calculate the mixing length (for dissipation)
[97]	111	dvar_dz = atmos_ocean_sign * &
	112	( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
	113	IF ( dvar_dz > 0.0 ) THEN
[57]	114	l_stable = 0.76 * SQRT( e(k,j,i) ) / &
[97]	115	SQRT( g / var_reference * dvar_dz ) + 1E-5
[57]	116	ELSE
[65]	117	l_stable = l_grid(k)
[57]	118	ENDIF
[1]	119	!
[65]	120	!-- Adjustment of the mixing length
	121	IF ( wall_adjustment ) THEN
[94]	122	l(k,j) = MIN( wall_adjustment_factor * &
	123	( zu(k) - zw(nzb_s_inner(j,i)) ), &
	124	l_grid(k), l_stable )
	125	ll(k,j) = MIN( wall_adjustment_factor * &
	126	( zu(k) - zw(nzb_s_inner(j,i)) ), &
	127	l_grid(k) )
[65]	128	ELSE
	129	l(k,j) = MIN( l_grid(k), l_stable )
	130	ll(k,j) = l_grid(k)
	131	ENDIF
	132	IF ( adjust_mixing_length .AND. prandtl_layer ) THEN
[94]	133	l(k,j) = MIN( l(k,j), kappa * &
	134	( zu(k) - zw(nzb_s_inner(j,i)) ) &
	135	/ phi_m )
	136	ll(k,j) = MIN( ll(k,j), kappa * &
	137	( zu(k) - zw(nzb_s_inner(j,i)) ) &
	138	/ phi_m )
[65]	139	ENDIF
[1]	140
[65]	141	ENDDO
[1]	142	ENDDO
[65]	143
[1]	144	!
[65]	145	!-- Calculate the tendency terms
	146	DO j = nys, nyn
	147	DO k = nzb_s_inner(j,i)+1, nzt
[1]	148
[65]	149	dissipation(k,j) = ( 0.19 + 0.74 * l(k,j) / ll(k,j) ) * &
	150	e(k,j,i) * SQRT( e(k,j,i) ) / l(k,j)
[1]	151
[65]	152	tend(k,j,i) = tend(k,j,i) &
[1]	153	+ ( &
	154	( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) ) &
	155	- ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) ) &
	156	) * ddx2 &
	157	+ ( &
	158	( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) ) &
	159	- ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) ) &
	160	) * ddy2 &
	161	+ ( &
	162	( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
	163	- ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k) &
	164	) * ddzw(k) &
	165	- dissipation(k,j)
	166
[65]	167	ENDDO
[1]	168	ENDDO
[65]	169
	170	!
	171	!-- Store dissipation if needed for calculating the sgs particle
	172	!-- velocities
[825]	173	IF ( use_sgs_for_particles .OR. wang_kernel ) THEN
[65]	174	DO j = nys, nyn
	175	DO k = nzb_s_inner(j,i)+1, nzt
	176	diss(k,j,i) = dissipation(k,j)
	177	ENDDO
	178	ENDDO
	179	ENDIF
	180
[1]	181	ENDDO
	182
[65]	183	ELSE
	184
	185	DO i = nxl, nxr
	186	DO j = nys, nyn
[1]	187	!
[65]	188	!-- First, calculate phi-function for eventually adjusting the &
	189	!-- mixing length to the prandtl mixing length
	190	IF ( adjust_mixing_length .AND. prandtl_layer ) THEN
	191	IF ( rif(j,i) >= 0.0 ) THEN
	192	phi_m = 1.0 + 5.0 * rif(j,i)
	193	ELSE
	194	phi_m = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rif(j,i) ) )
	195	ENDIF
	196	ENDIF
	197
	198	DO k = nzb_s_inner(j,i)+1, nzt
	199	!
	200	!-- Calculate the mixing length (for dissipation)
[97]	201	dvar_dz = atmos_ocean_sign * &
	202	( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
	203	IF ( dvar_dz > 0.0 ) THEN
[65]	204	l_stable = 0.76 * SQRT( e(k,j,i) ) / &
[97]	205	SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5
[65]	206	ELSE
	207	l_stable = l_grid(k)
	208	ENDIF
	209	!
	210	!-- Adjustment of the mixing length
	211	IF ( wall_adjustment ) THEN
[94]	212	l(k,j) = MIN( wall_adjustment_factor * &
	213	( zu(k) - zw(nzb_s_inner(j,i)) ), &
	214	l_grid(k), l_stable )
	215	ll(k,j) = MIN( wall_adjustment_factor * &
	216	( zu(k) - zw(nzb_s_inner(j,i)) ), &
	217	l_grid(k) )
[65]	218	ELSE
	219	l(k,j) = MIN( l_grid(k), l_stable )
	220	ll(k,j) = l_grid(k)
	221	ENDIF
	222	IF ( adjust_mixing_length .AND. prandtl_layer ) THEN
[94]	223	l(k,j) = MIN( l(k,j), kappa * &
	224	( zu(k) - zw(nzb_s_inner(j,i)) ) &
	225	/ phi_m )
	226	ll(k,j) = MIN( ll(k,j), kappa * &
	227	( zu(k) - zw(nzb_s_inner(j,i)) ) &
	228	/ phi_m )
[65]	229	ENDIF
	230
	231	ENDDO
	232	ENDDO
	233
	234	!
	235	!-- Calculate the tendency terms
[1]	236	DO j = nys, nyn
[19]	237	DO k = nzb_s_inner(j,i)+1, nzt
[65]	238
	239	dissipation(k,j) = ( 0.19 + 0.74 * l(k,j) / ll(k,j) ) * &
	240	e(k,j,i) * SQRT( e(k,j,i) ) / l(k,j)
	241
	242	tend(k,j,i) = tend(k,j,i) &
	243	+ ( &
	244	( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) ) &
	245	- ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) ) &
	246	) * ddx2 &
	247	+ ( &
	248	( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) ) &
	249	- ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) ) &
	250	) * ddy2 &
	251	+ ( &
	252	( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
	253	- ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k) &
	254	) * ddzw(k) &
	255	- dissipation(k,j)
	256
[1]	257	ENDDO
	258	ENDDO
	259
[65]	260	!
	261	!-- Store dissipation if needed for calculating the sgs particle
	262	!-- velocities
[825]	263	IF ( use_sgs_for_particles .OR. wang_kernel ) THEN
[65]	264	DO j = nys, nyn
	265	DO k = nzb_s_inner(j,i)+1, nzt
	266	diss(k,j,i) = dissipation(k,j)
	267	ENDDO
	268	ENDDO
	269	ENDIF
[1]	270
[65]	271	ENDDO
	272
	273	ENDIF
	274
[1]	275	!
	276	!-- Boundary condition for dissipation
[825]	277	IF ( use_sgs_for_particles .OR. wang_kernel ) THEN
[1]	278	DO i = nxl, nxr
	279	DO j = nys, nyn
	280	diss(nzb_s_inner(j,i),j,i) = diss(nzb_s_inner(j,i)+1,j,i)
	281	ENDDO
	282	ENDDO
	283	ENDIF
	284
	285	END SUBROUTINE diffusion_e
	286
	287
	288	!------------------------------------------------------------------------------!
	289	! Call for grid point i,j
	290	!------------------------------------------------------------------------------!
[1001]	291	SUBROUTINE diffusion_e_ij( i, j, var, var_reference )
[1]	292
[1001]	293	USE arrays_3d
[1]	294	USE control_parameters
	295	USE grid_variables
	296	USE indices
	297	USE particle_attributes
	298
	299	IMPLICIT NONE
	300
[1001]	301	INTEGER :: i, j, k
	302	REAL :: dvar_dz, l_stable, phi_m, var_reference
[1]	303
[1010]	304	#if defined( __nopointer )
	305	REAL, DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: var
	306	#else
[1001]	307	REAL, DIMENSION(:,:,:), POINTER :: var
[1010]	308	#endif
[1001]	309	REAL, DIMENSION(nzb+1:nzt) :: dissipation, l, ll
[1]	310
[1001]	311
[1]	312	!
	313	!-- First, calculate phi-function for eventually adjusting the mixing length
	314	!-- to the prandtl mixing length
	315	IF ( adjust_mixing_length .AND. prandtl_layer ) THEN
	316	IF ( rif(j,i) >= 0.0 ) THEN
	317	phi_m = 1.0 + 5.0 * rif(j,i)
	318	ELSE
	319	phi_m = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rif(j,i) ) )
	320	ENDIF
	321	ENDIF
	322
	323	!
	324	!-- Calculate the mixing length (for dissipation)
[19]	325	DO k = nzb_s_inner(j,i)+1, nzt
[97]	326	dvar_dz = atmos_ocean_sign * &
	327	( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
	328	IF ( dvar_dz > 0.0 ) THEN
	329	IF ( use_reference ) THEN
[57]	330	l_stable = 0.76 * SQRT( e(k,j,i) ) / &
[97]	331	SQRT( g / var_reference * dvar_dz ) + 1E-5
[57]	332	ELSE
	333	l_stable = 0.76 * SQRT( e(k,j,i) ) / &
[97]	334	SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5
[57]	335	ENDIF
[1]	336	ELSE
	337	l_stable = l_grid(k)
	338	ENDIF
	339	!
	340	!-- Adjustment of the mixing length
	341	IF ( wall_adjustment ) THEN
[94]	342	l(k) = MIN( wall_adjustment_factor * &
	343	( zu(k) - zw(nzb_s_inner(j,i)) ), l_grid(k), &
	344	l_stable )
	345	ll(k) = MIN( wall_adjustment_factor * &
	346	( zu(k) - zw(nzb_s_inner(j,i)) ), l_grid(k) )
[1]	347	ELSE
	348	l(k) = MIN( l_grid(k), l_stable )
	349	ll(k) = l_grid(k)
	350	ENDIF
	351	IF ( adjust_mixing_length .AND. prandtl_layer ) THEN
[94]	352	l(k) = MIN( l(k), kappa * &
	353	( zu(k) - zw(nzb_s_inner(j,i)) ) / phi_m )
	354	ll(k) = MIN( ll(k), kappa * &
	355	( zu(k) - zw(nzb_s_inner(j,i)) ) / phi_m )
[1]	356	ENDIF
	357
	358	!
	359	!-- Calculate the tendency term
	360	dissipation(k) = ( 0.19 + 0.74 * l(k) / ll(k) ) * e(k,j,i) * &
	361	SQRT( e(k,j,i) ) / l(k)
	362
	363	tend(k,j,i) = tend(k,j,i) &
	364	+ ( &
	365	( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) ) &
	366	- ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) ) &
	367	) * ddx2 &
	368	+ ( &
	369	( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) ) &
	370	- ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) ) &
	371	) * ddy2 &
	372	+ ( &
	373	( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
	374	- ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k) &
	375	) * ddzw(k) &
	376	- dissipation(k)
	377
	378	ENDDO
	379
	380	!
	381	!-- Store dissipation if needed for calculating the sgs particle velocities
[825]	382	IF ( use_sgs_for_particles .OR. wang_kernel ) THEN
[19]	383	DO k = nzb_s_inner(j,i)+1, nzt
[1]	384	diss(k,j,i) = dissipation(k)
	385	ENDDO
	386	!
	387	!-- Boundary condition for dissipation
	388	diss(nzb_s_inner(j,i),j,i) = diss(nzb_s_inner(j,i)+1,j,i)
	389	ENDIF
	390
	391	END SUBROUTINE diffusion_e_ij
	392
	393	END MODULE diffusion_e_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |