Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |