Home

Context Navigation

diffusion_v.f90 @ 219

Last change on this file since 219 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.6 KB

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

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

Context Navigation

source: palm/trunk/SOURCE/diffusion_v.f90 @ 219

Download in other formats: