Home

Context Navigation

diffusion_u.f90 @ 250

Last change on this file since 250 was 110, checked in by raasch, 17 years ago

New:
---
Allows runs for a coupled atmosphere-ocean LES,
coupling frequency is controlled by new d3par-parameter dt_coupling,
the coupling mode (atmosphere_to_ocean or ocean_to_atmosphere) for the
respective processes is read from environment variable coupling_mode,
which is set by the mpiexec-command,
communication between the two models is done using the intercommunicator
comm_inter,
local files opened by the ocean model get the additional suffic "_O".
Assume saturation at k=nzb_s_inner(j,i) for atmosphere coupled to ocean.

A momentum flux can be set as top boundary condition using the new
inipar parameter top_momentumflux_u|v.

Non-cyclic boundary conditions can be used along all horizontal directions.

Quantities w*p* and w"e can be output as vertical profiles.

Initial profiles are reset to constant profiles in case that initializing_actions /= 'set_constant_profiles'. (init_rankine)

Optionally calculate km and kh from initial TKE e_init.

Changed:

Remaining variables iran changed to iran_part (advec_particles, init_particles).

In case that the presure solver is not called for every Runge-Kutta substep
(call_psolver_at_all_substeps = .F.), it is called after the first substep
instead of the last. In that case, random perturbations are also added to the
velocity field after the first substep.

Initialization of km,kh = 0.00001 for ocean = .T. (for ocean = .F. it remains 0.01).

Allow data_output_pr= q, wq, w"q", w*q* for humidity = .T. (instead of cloud_physics = .T.).

Errors:

Bugs from code parts for non-cyclic boundary conditions are removed: loops for
u and v are starting from index nxlu, nysv, respectively. The radiation boundary
condition is used for every Runge-Kutta substep. Velocity phase speeds for
the radiation boundary conditions are calculated for the first Runge-Kutta
substep only and reused for the further substeps. New arrays c_u, c_v, and c_w
are defined for this purpose. Several index errors are removed from the
radiation boundary condition code parts. Upper bounds for calculating
u_0 and v_0 (in production_e) are nxr+1 and nyn+1 because otherwise these
values are not available in case of non-cyclic boundary conditions.

+dots_num_palm in module user, +module netcdf_control in user_init (both in user_interface)

Bugfix: wrong sign removed from the buoyancy production term in the case use_reference = .T. (production_e)

Bugfix: Error message concerning output of particle concentration (pc) modified (check_parameters).

Bugfix: Rayleigh damping for ocean fixed.

Property svn:keywords set to Id

File size: 17.7 KB

Rev	Line
[1]	1	MODULE diffusion_u_mod
	2
	3	!------------------------------------------------------------------------------!
	4	! Actual revisions:
	5	! -----------------
[110]	6	!
[106]	7	!
[1]	8	! Former revisions:
	9	! -----------------
[3]	10	! $Id: diffusion_u.f90 110 2007-10-05 05:13:14Z letzel $
[39]	11	!
[110]	12	! 106 2007-08-16 14:30:26Z raasch
	13	! Momentumflux at top (uswst) included as boundary condition,
	14	! i loop is starting from nxlu (needed for non-cyclic boundary conditions)
	15	!
[77]	16	! 75 2007-03-22 09:54:05Z raasch
	17	! Wall functions now include diabatic conditions, call of routine wall_fluxes,
	18	! z0 removed from argument list, uxrp eliminated
	19	!
[39]	20	! 20 2007-02-26 00:12:32Z raasch
	21	! Bugfix: ddzw dimensioned 1:nzt"+1"
	22	!
[3]	23	! RCS Log replace by Id keyword, revision history cleaned up
	24	!
[1]	25	! Revision 1.15 2006/02/23 10:35:35 raasch
	26	! nzb_2d replaced by nzb_u_outer in horizontal diffusion and by nzb_u_inner
	27	! or nzb_diff_u, respectively, in vertical diffusion,
	28	! wall functions added for north and south walls, +z0 in argument list,
	29	! terms containing w(k-1,..) are removed from the Prandtl-layer equation
	30	! because they cause errors at the edges of topography
	31	! WARNING: loops containing the MAX function are still not properly vectorized!
	32	!
	33	! Revision 1.1 1997/09/12 06:23:51 raasch
	34	! Initial revision
	35	!
	36	!
	37	! Description:
	38	! ------------
	39	! Diffusion term of the u-component
[51]	40	! To do: additional damping (needed for non-cyclic bc) causes bad vectorization
	41	! and slows down the speed on NEC about 5-10%
[1]	42	!------------------------------------------------------------------------------!
	43
[56]	44	USE wall_fluxes_mod
	45
[1]	46	PRIVATE
	47	PUBLIC diffusion_u
	48
	49	INTERFACE diffusion_u
	50	MODULE PROCEDURE diffusion_u
	51	MODULE PROCEDURE diffusion_u_ij
	52	END INTERFACE diffusion_u
	53
	54	CONTAINS
	55
	56
	57	!------------------------------------------------------------------------------!
	58	! Call for all grid points
	59	!------------------------------------------------------------------------------!
[102]	60	SUBROUTINE diffusion_u( ddzu, ddzw, km, km_damp_y, tend, u, usws, uswst, &
	61	v, w )
[1]	62
	63	USE control_parameters
	64	USE grid_variables
	65	USE indices
	66
	67	IMPLICIT NONE
	68
	69	INTEGER :: i, j, k
[51]	70	REAL :: kmym_x, kmym_y, kmyp_x, kmyp_y, kmzm, kmzp
[20]	71	REAL :: ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_y(nys-1:nyn+1)
[1]	72	REAL :: tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
[102]	73	REAL, DIMENSION(:,:), POINTER :: usws, uswst
[1]	74	REAL, DIMENSION(:,:,:), POINTER :: km, u, v, w
[75]	75	REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: usvs
[1]	76
[56]	77	!
	78	!-- First calculate horizontal momentum flux u'v' at vertical walls,
	79	!-- if neccessary
	80	IF ( topography /= 'flat' ) THEN
[75]	81	CALL wall_fluxes( usvs, 1.0, 0.0, 0.0, 0.0, nzb_u_inner, &
[56]	82	nzb_u_outer, wall_u )
	83	ENDIF
	84
[106]	85	DO i = nxlu, nxr
[1]	86	DO j = nys,nyn
	87	!
	88	!-- Compute horizontal diffusion
	89	DO k = nzb_u_outer(j,i)+1, nzt
	90	!
	91	!-- Interpolate eddy diffusivities on staggered gridpoints
	92	kmyp_x = 0.25 * &
	93	( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
	94	kmym_x = 0.25 * &
	95	( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
	96	kmyp_y = kmyp_x
	97	kmym_y = kmym_x
	98	!
	99	!-- Increase diffusion at the outflow boundary in case of
	100	!-- non-cyclic lateral boundaries. Damping is only needed for
	101	!-- velocity components parallel to the outflow boundary in
	102	!-- the direction normal to the outflow boundary.
	103	IF ( bc_ns /= 'cyclic' ) THEN
	104	kmyp_y = MAX( kmyp_y, km_damp_y(j) )
	105	kmym_y = MAX( kmym_y, km_damp_y(j) )
	106	ENDIF
	107
	108	tend(k,j,i) = tend(k,j,i) &
	109	& + 2.0 * ( &
	110	& km(k,j,i) * ( u(k,j,i+1) - u(k,j,i) ) &
	111	& - km(k,j,i-1) * ( u(k,j,i) - u(k,j,i-1) ) &
	112	& ) * ddx2 &
	113	& + ( kmyp_y * ( u(k,j+1,i) - u(k,j,i) ) * ddy &
	114	& + kmyp_x * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx &
	115	& - kmym_y * ( u(k,j,i) - u(k,j-1,i) ) * ddy &
	116	& - kmym_x * ( v(k,j,i) - v(k,j,i-1) ) * ddx &
	117	& ) * ddy
	118	ENDDO
	119
	120	!
	121	!-- Wall functions at the north and south walls, respectively
	122	IF ( wall_u(j,i) /= 0.0 ) THEN
[51]	123
[1]	124	DO k = nzb_u_inner(j,i)+1, nzb_u_outer(j,i)
	125	kmyp_x = 0.25 * &
	126	( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
	127	kmym_x = 0.25 * &
	128	( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
	129	kmyp_y = kmyp_x
	130	kmym_y = kmym_x
	131	!
	132	!-- Increase diffusion at the outflow boundary in case of
	133	!-- non-cyclic lateral boundaries. Damping is only needed for
	134	!-- velocity components parallel to the outflow boundary in
	135	!-- the direction normal to the outflow boundary.
	136	IF ( bc_ns /= 'cyclic' ) THEN
	137	kmyp_y = MAX( kmyp_y, km_damp_y(j) )
	138	kmym_y = MAX( kmym_y, km_damp_y(j) )
	139	ENDIF
	140
	141	tend(k,j,i) = tend(k,j,i) &
	142	+ 2.0 * ( &
	143	km(k,j,i) * ( u(k,j,i+1) - u(k,j,i) ) &
	144	- km(k,j,i-1) * ( u(k,j,i) - u(k,j,i-1) ) &
	145	) * ddx2 &
	146	+ ( fyp(j,i) * ( &
	147	kmyp_y * ( u(k,j+1,i) - u(k,j,i) ) * ddy &
	148	+ kmyp_x * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx &
	149	) &
	150	- fym(j,i) * ( &
	151	kmym_y * ( u(k,j,i) - u(k,j-1,i) ) * ddy &
	152	+ kmym_x * ( v(k,j,i) - v(k,j,i-1) ) * ddx &
	153	) &
[56]	154	+ wall_u(j,i) * usvs(k,j,i) &
[1]	155	) * ddy
	156	ENDDO
	157	ENDIF
	158
	159	!
	160	!-- Compute vertical diffusion. In case of simulating a Prandtl layer,
	161	!-- index k starts at nzb_u_inner+2.
[102]	162	DO k = nzb_diff_u(j,i), nzt_diff
[1]	163	!
	164	!-- Interpolate eddy diffusivities on staggered gridpoints
	165	kmzp = 0.25 * &
	166	( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
	167	kmzm = 0.25 * &
	168	( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
	169
	170	tend(k,j,i) = tend(k,j,i) &
	171	& + ( kmzp * ( ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1) &
	172	& + ( w(k,j,i) - w(k,j,i-1) ) * ddx &
	173	& ) &
	174	& - kmzm * ( ( u(k,j,i) - u(k-1,j,i) ) * ddzu(k) &
	175	& + ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx &
	176	& ) &
	177	& ) * ddzw(k)
	178	ENDDO
	179
	180	!
	181	!-- Vertical diffusion at the first grid point above the surface,
	182	!-- if the momentum flux at the bottom is given by the Prandtl law or
	183	!-- if it is prescribed by the user.
	184	!-- Difference quotient of the momentum flux is not formed over half
	185	!-- of the grid spacing (2.0*ddzw(k)) any more, since the comparison
	186	!-- with other (LES) modell showed that the values of the momentum
	187	!-- flux becomes too large in this case.
	188	!-- The term containing w(k-1,..) (see above equation) is removed here
	189	!-- because the vertical velocity is assumed to be zero at the surface.
	190	IF ( use_surface_fluxes ) THEN
	191	k = nzb_u_inner(j,i)+1
	192	!
	193	!-- Interpolate eddy diffusivities on staggered gridpoints
	194	kmzp = 0.25 * &
	195	( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
	196	kmzm = 0.25 * &
	197	( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
	198
	199	tend(k,j,i) = tend(k,j,i) &
	200	& + ( kmzp * ( w(k,j,i) - w(k,j,i-1) ) * ddx &
	201	& ) * ddzw(k) &
[102]	202	& + ( kmzp * ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1) &
[1]	203	& + usws(j,i) &
	204	& ) * ddzw(k)
	205	ENDIF
	206
[102]	207	!
	208	!-- Vertical diffusion at the first gridpoint below the top boundary,
	209	!-- if the momentum flux at the top is prescribed by the user
[103]	210	IF ( use_top_fluxes .AND. constant_top_momentumflux ) THEN
[102]	211	k = nzt
	212	!
	213	!-- Interpolate eddy diffusivities on staggered gridpoints
	214	kmzp = 0.25 * &
	215	( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
	216	kmzm = 0.25 * &
	217	( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
	218
	219	tend(k,j,i) = tend(k,j,i) &
	220	& - ( kmzm * ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx &
	221	& ) * ddzw(k) &
	222	& + ( -uswst(j,i) &
	223	& - kmzm * ( u(k,j,i) - u(k-1,j,i) ) * ddzu(k) &
	224	& ) * ddzw(k)
	225	ENDIF
	226
[1]	227	ENDDO
	228	ENDDO
	229
	230	END SUBROUTINE diffusion_u
	231
	232
	233	!------------------------------------------------------------------------------!
	234	! Call for grid point i,j
	235	!------------------------------------------------------------------------------!
	236	SUBROUTINE diffusion_u_ij( i, j, ddzu, ddzw, km, km_damp_y, tend, u, usws, &
[102]	237	uswst, v, w )
[1]	238
	239	USE control_parameters
	240	USE grid_variables
	241	USE indices
	242
	243	IMPLICIT NONE
	244
	245	INTEGER :: i, j, k
[51]	246	REAL :: kmym_x, kmym_y, kmyp_x, kmyp_y, kmzm, kmzp
[20]	247	REAL :: ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_y(nys-1:nyn+1)
[1]	248	REAL :: tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
[51]	249	REAL, DIMENSION(nzb:nzt+1) :: usvs
[102]	250	REAL, DIMENSION(:,:), POINTER :: usws, uswst
[1]	251	REAL, DIMENSION(:,:,:), POINTER :: km, u, v, w
	252
	253	!
	254	!-- Compute horizontal diffusion
	255	DO k = nzb_u_outer(j,i)+1, nzt
	256	!
	257	!-- Interpolate eddy diffusivities on staggered gridpoints
	258	kmyp_x = 0.25 * ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
	259	kmym_x = 0.25 * ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
	260	kmyp_y = kmyp_x
	261	kmym_y = kmym_x
	262
	263	!
	264	!-- Increase diffusion at the outflow boundary in case of non-cyclic
	265	!-- lateral boundaries. Damping is only needed for velocity components
	266	!-- parallel to the outflow boundary in the direction normal to the
	267	!-- outflow boundary.
	268	IF ( bc_ns /= 'cyclic' ) THEN
	269	kmyp_y = MAX( kmyp_y, km_damp_y(j) )
	270	kmym_y = MAX( kmym_y, km_damp_y(j) )
	271	ENDIF
	272
	273	tend(k,j,i) = tend(k,j,i) &
	274	& + 2.0 * ( &
	275	& km(k,j,i) * ( u(k,j,i+1) - u(k,j,i) ) &
	276	& - km(k,j,i-1) * ( u(k,j,i) - u(k,j,i-1) ) &
	277	& ) * ddx2 &
	278	& + ( kmyp_y * ( u(k,j+1,i) - u(k,j,i) ) * ddy &
	279	& + kmyp_x * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx &
	280	& - kmym_y * ( u(k,j,i) - u(k,j-1,i) ) * ddy &
	281	& - kmym_x * ( v(k,j,i) - v(k,j,i-1) ) * ddx &
	282	& ) * ddy
	283	ENDDO
	284
	285	!
	286	!-- Wall functions at the north and south walls, respectively
	287	IF ( wall_u(j,i) .NE. 0.0 ) THEN
[51]	288
	289	!
	290	!-- Calculate the horizontal momentum flux u'v'
	291	CALL wall_fluxes( i, j, nzb_u_inner(j,i)+1, nzb_u_outer(j,i), &
	292	usvs, 1.0, 0.0, 0.0, 0.0 )
	293
[1]	294	DO k = nzb_u_inner(j,i)+1, nzb_u_outer(j,i)
	295	kmyp_x = 0.25 * ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
	296	kmym_x = 0.25 * ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
	297	kmyp_y = kmyp_x
	298	kmym_y = kmym_x
	299	!
	300	!-- Increase diffusion at the outflow boundary in case of
	301	!-- non-cyclic lateral boundaries. Damping is only needed for
	302	!-- velocity components parallel to the outflow boundary in
	303	!-- the direction normal to the outflow boundary.
	304	IF ( bc_ns /= 'cyclic' ) THEN
	305	kmyp_y = MAX( kmyp_y, km_damp_y(j) )
	306	kmym_y = MAX( kmym_y, km_damp_y(j) )
	307	ENDIF
	308
	309	tend(k,j,i) = tend(k,j,i) &
	310	+ 2.0 * ( &
	311	km(k,j,i) * ( u(k,j,i+1) - u(k,j,i) ) &
	312	- km(k,j,i-1) * ( u(k,j,i) - u(k,j,i-1) ) &
	313	) * ddx2 &
	314	+ ( fyp(j,i) * ( &
	315	kmyp_y * ( u(k,j+1,i) - u(k,j,i) ) * ddy &
	316	+ kmyp_x * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx &
	317	) &
	318	- fym(j,i) * ( &
	319	kmym_y * ( u(k,j,i) - u(k,j-1,i) ) * ddy &
	320	+ kmym_x * ( v(k,j,i) - v(k,j,i-1) ) * ddx &
	321	) &
[51]	322	+ wall_u(j,i) * usvs(k) &
[1]	323	) * ddy
	324	ENDDO
	325	ENDIF
	326
	327	!
	328	!-- Compute vertical diffusion. In case of simulating a Prandtl layer,
	329	!-- index k starts at nzb_u_inner+2.
[102]	330	DO k = nzb_diff_u(j,i), nzt_diff
[1]	331	!
	332	!-- Interpolate eddy diffusivities on staggered gridpoints
	333	kmzp = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
	334	kmzm = 0.25 * ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
	335
	336	tend(k,j,i) = tend(k,j,i) &
	337	& + ( kmzp * ( ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1) &
	338	& + ( w(k,j,i) - w(k,j,i-1) ) * ddx &
	339	& ) &
	340	& - kmzm * ( ( u(k,j,i) - u(k-1,j,i) ) * ddzu(k) &
	341	& + ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx &
	342	& ) &
	343	& ) * ddzw(k)
	344	ENDDO
	345
	346	!
	347	!-- Vertical diffusion at the first grid point above the surface, if the
	348	!-- momentum flux at the bottom is given by the Prandtl law or if it is
	349	!-- prescribed by the user.
	350	!-- Difference quotient of the momentum flux is not formed over half of
	351	!-- the grid spacing (2.0*ddzw(k)) any more, since the comparison with
	352	!-- other (LES) modell showed that the values of the momentum flux becomes
	353	!-- too large in this case.
	354	!-- The term containing w(k-1,..) (see above equation) is removed here
	355	!-- because the vertical velocity is assumed to be zero at the surface.
	356	IF ( use_surface_fluxes ) THEN
	357	k = nzb_u_inner(j,i)+1
	358	!
	359	!-- Interpolate eddy diffusivities on staggered gridpoints
	360	kmzp = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
	361	kmzm = 0.25 * ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
	362
	363	tend(k,j,i) = tend(k,j,i) &
	364	& + ( kmzp * ( w(k,j,i) - w(k,j,i-1) ) * ddx &
	365	& ) * ddzw(k) &
[102]	366	& + ( kmzp * ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1) &
[1]	367	& + usws(j,i) &
	368	& ) * ddzw(k)
	369	ENDIF
	370
[102]	371	!
	372	!-- Vertical diffusion at the first gridpoint below the top boundary,
	373	!-- if the momentum flux at the top is prescribed by the user
[103]	374	IF ( use_top_fluxes .AND. constant_top_momentumflux ) THEN
[102]	375	k = nzt
	376	!
	377	!-- Interpolate eddy diffusivities on staggered gridpoints
	378	kmzp = 0.25 * &
	379	( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
	380	kmzm = 0.25 * &
	381	( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
	382
	383	tend(k,j,i) = tend(k,j,i) &
	384	& - ( kmzm * ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx &
	385	& ) * ddzw(k) &
	386	& + ( -uswst(j,i) &
	387	& - kmzm * ( u(k,j,i) - u(k-1,j,i) ) * ddzu(k) &
	388	& ) * ddzw(k)
	389	ENDIF
	390
[1]	391	END SUBROUTINE diffusion_u_ij
	392
	393	END MODULE diffusion_u_mod

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

Context Navigation

source: palm/trunk/SOURCE/diffusion_u.f90 @ 250

Download in other formats: