Home

Context Navigation

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

Last change on this file since 2045 was 2038, checked in by knoop, 8 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 27.1 KB

Rev	Line
[1873]	1	!> @file diffusion_v.f90
[2000]	2	!------------------------------------------------------------------------------!
[1036]	3	! This file is part of PALM.
	4	!
[2000]	5	! PALM is free software: you can redistribute it and/or modify it under the
	6	! terms of the GNU General Public License as published by the Free Software
	7	! Foundation, either version 3 of the License, or (at your option) any later
	8	! version.
[1036]	9	!
	10	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
	11	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
	12	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
	13	!
	14	! You should have received a copy of the GNU General Public License along with
	15	! PALM. If not, see <http://www.gnu.org/licenses/>.
	16	!
[1818]	17	! Copyright 1997-2016 Leibniz Universitaet Hannover
[2000]	18	!------------------------------------------------------------------------------!
[1036]	19	!
[484]	20	! Current revisions:
[1]	21	! -----------------
[1341]	22	!
[2038]	23	!
[1321]	24	! Former revisions:
	25	! -----------------
	26	! $Id: diffusion_v.f90 2038 2016-10-26 11:16:56Z knoop $
	27	!
[2038]	28	! 2037 2016-10-26 11:15:40Z knoop
	29	! Anelastic approximation implemented
	30	!
[2001]	31	! 2000 2016-08-20 18:09:15Z knoop
	32	! Forced header and separation lines into 80 columns
	33	!
[1874]	34	! 1873 2016-04-18 14:50:06Z maronga
	35	! Module renamed (removed _mod)
	36	!
	37	!
[1851]	38	! 1850 2016-04-08 13:29:27Z maronga
	39	! Module renamed
	40	!
	41	!
[1741]	42	! 1740 2016-01-13 08:19:40Z raasch
	43	! unnecessary calculations of kmzm and kmzp in wall bounded parts removed
	44	!
[1683]	45	! 1682 2015-10-07 23:56:08Z knoop
	46	! Code annotations made doxygen readable
	47	!
[1341]	48	! 1340 2014-03-25 19:45:13Z kanani
	49	! REAL constants defined as wp-kind
	50	!
[1321]	51	! 1320 2014-03-20 08:40:49Z raasch
[1320]	52	! ONLY-attribute added to USE-statements,
	53	! kind-parameters added to all INTEGER and REAL declaration statements,
	54	! kinds are defined in new module kinds,
	55	! revision history before 2012 removed,
	56	! comment fields (!:) to be used for variable explanations added to
	57	! all variable declaration statements
[1321]	58	!
[1258]	59	! 1257 2013-11-08 15:18:40Z raasch
	60	! openacc loop and loop vector clauses removed, declare create moved after
	61	! the FORTRAN declaration statement
	62	!
[1132]	63	! 1128 2013-04-12 06:19:32Z raasch
	64	! loop index bounds in accelerator version replaced by i_left, i_right, j_south,
	65	! j_north
	66	!
[1037]	67	! 1036 2012-10-22 13:43:42Z raasch
	68	! code put under GPL (PALM 3.9)
	69	!
[1017]	70	! 1015 2012-09-27 09:23:24Z raasch
	71	! accelerator version (*_acc) added
	72	!
[1002]	73	! 1001 2012-09-13 14:08:46Z raasch
	74	! arrays comunicated by module instead of parameter list
	75	!
[979]	76	! 978 2012-08-09 08:28:32Z fricke
	77	! outflow damping layer removed
	78	! kmxm_x/_y and kmxp_x/_y change to kmxm and kmxp
	79	!
[1]	80	! Revision 1.1 1997/09/12 06:24:01 raasch
	81	! Initial revision
	82	!
	83	!
	84	! Description:
	85	! ------------
[1682]	86	!> Diffusion term of the v-component
[1]	87	!------------------------------------------------------------------------------!
[1682]	88	MODULE diffusion_v_mod
	89
[1]	90
[56]	91	USE wall_fluxes_mod
	92
[1]	93	PRIVATE
[1015]	94	PUBLIC diffusion_v, diffusion_v_acc
[1]	95
	96	INTERFACE diffusion_v
	97	MODULE PROCEDURE diffusion_v
	98	MODULE PROCEDURE diffusion_v_ij
	99	END INTERFACE diffusion_v
	100
[1015]	101	INTERFACE diffusion_v_acc
	102	MODULE PROCEDURE diffusion_v_acc
	103	END INTERFACE diffusion_v_acc
	104
[1]	105	CONTAINS
	106
	107
	108	!------------------------------------------------------------------------------!
[1682]	109	! Description:
	110	! ------------
	111	!> Call for all grid points
[1]	112	!------------------------------------------------------------------------------!
[1001]	113	SUBROUTINE diffusion_v
[1]	114
[1320]	115	USE arrays_3d, &
[2037]	116	ONLY: ddzu, ddzw, km, tend, u, v, vsws, vswst, w, &
	117	drho_air, rho_air_zw
[1320]	118
	119	USE control_parameters, &
	120	ONLY: constant_top_momentumflux, topography, use_surface_fluxes, &
	121	use_top_fluxes
	122
	123	USE grid_variables, &
	124	ONLY: ddx, ddy, ddy2, fxm, fxp, wall_v
	125
	126	USE indices, &
	127	ONLY: nxl, nxr, nyn, nys, nysv, nzb, nzb_diff_v, nzb_v_inner, &
	128	nzb_v_outer, nzt, nzt_diff
	129
	130	USE kinds
[1]	131
	132	IMPLICIT NONE
	133
[1682]	134	INTEGER(iwp) :: i !<
	135	INTEGER(iwp) :: j !<
	136	INTEGER(iwp) :: k !<
	137	REAL(wp) :: kmxm !<
	138	REAL(wp) :: kmxp !<
	139	REAL(wp) :: kmzm !<
	140	REAL(wp) :: kmzp !<
[1001]	141
[1682]	142	REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: vsus !<
[1]	143
[56]	144	!
	145	!-- First calculate horizontal momentum flux v'u' at vertical walls,
	146	!-- if neccessary
	147	IF ( topography /= 'flat' ) THEN
[1320]	148	CALL wall_fluxes( vsus, 0.0_wp, 1.0_wp, 0.0_wp, 0.0_wp, nzb_v_inner, &
[56]	149	nzb_v_outer, wall_v )
	150	ENDIF
	151
[1]	152	DO i = nxl, nxr
[106]	153	DO j = nysv, nyn
[1]	154	!
	155	!-- Compute horizontal diffusion
	156	DO k = nzb_v_outer(j,i)+1, nzt
	157	!
	158	!-- Interpolate eddy diffusivities on staggered gridpoints
[1340]	159	kmxp = 0.25_wp * &
[978]	160	( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
[1340]	161	kmxm = 0.25_wp * &
[978]	162	( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
[1]	163
[1320]	164	tend(k,j,i) = tend(k,j,i) &
	165	& + ( kmxp * ( v(k,j,i+1) - v(k,j,i) ) * ddx &
	166	& + kmxp * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy &
	167	& - kmxm * ( v(k,j,i) - v(k,j,i-1) ) * ddx &
	168	& - kmxm * ( u(k,j,i) - u(k,j-1,i) ) * ddy &
	169	& ) * ddx &
[1340]	170	& + 2.0_wp * ( &
[1320]	171	& km(k,j,i) * ( v(k,j+1,i) - v(k,j,i) ) &
	172	& - km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) ) &
[1340]	173	& ) * ddy2
[1]	174	ENDDO
	175
	176	!
	177	!-- Wall functions at the left and right walls, respectively
[1340]	178	IF ( wall_v(j,i) /= 0.0_wp ) THEN
[51]	179
[1]	180	DO k = nzb_v_inner(j,i)+1, nzb_v_outer(j,i)
[1340]	181	kmxp = 0.25_wp * &
[978]	182	( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
[1340]	183	kmxm = 0.25_wp * &
[978]	184	( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
	185
[1]	186	tend(k,j,i) = tend(k,j,i) &
[1340]	187	+ 2.0_wp * ( &
[1]	188	km(k,j,i) * ( v(k,j+1,i) - v(k,j,i) ) &
	189	- km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) ) &
[1340]	190	) * ddy2 &
[1]	191	+ ( fxp(j,i) * ( &
[978]	192	kmxp * ( v(k,j,i+1) - v(k,j,i) ) * ddx &
	193	+ kmxp * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy &
[1]	194	) &
	195	- fxm(j,i) * ( &
[978]	196	kmxm * ( v(k,j,i) - v(k,j,i-1) ) * ddx &
	197	+ kmxm * ( u(k,j,i) - u(k,j-1,i) ) * ddy &
[1]	198	) &
[56]	199	+ wall_v(j,i) * vsus(k,j,i) &
[1]	200	) * ddx
	201	ENDDO
	202	ENDIF
	203
	204	!
	205	!-- Compute vertical diffusion. In case of simulating a Prandtl
	206	!-- layer, index k starts at nzb_v_inner+2.
[102]	207	DO k = nzb_diff_v(j,i), nzt_diff
[1]	208	!
	209	!-- Interpolate eddy diffusivities on staggered gridpoints
[1340]	210	kmzp = 0.25_wp * &
[1]	211	( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
[1340]	212	kmzm = 0.25_wp * &
[1]	213	( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
	214
[1320]	215	tend(k,j,i) = tend(k,j,i) &
	216	& + ( kmzp * ( ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1) &
	217	& + ( w(k,j,i) - w(k,j-1,i) ) * ddy &
[2037]	218	& ) * rho_air_zw(k) &
[1320]	219	& - kmzm * ( ( v(k,j,i) - v(k-1,j,i) ) * ddzu(k) &
	220	& + ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy &
[2037]	221	& ) * rho_air_zw(k-1) &
	222	& ) * ddzw(k) * drho_air(k)
[1]	223	ENDDO
	224
	225	!
	226	!-- Vertical diffusion at the first grid point above the surface,
	227	!-- if the momentum flux at the bottom is given by the Prandtl law
	228	!-- or if it is prescribed by the user.
	229	!-- Difference quotient of the momentum flux is not formed over
	230	!-- half of the grid spacing (2.0*ddzw(k)) any more, since the
[1320]	231	!-- comparison with other (LES) models showed that the values of
[1]	232	!-- the momentum flux becomes too large in this case.
	233	!-- The term containing w(k-1,..) (see above equation) is removed here
	234	!-- because the vertical velocity is assumed to be zero at the surface.
	235	IF ( use_surface_fluxes ) THEN
	236	k = nzb_v_inner(j,i)+1
	237	!
	238	!-- Interpolate eddy diffusivities on staggered gridpoints
[1340]	239	kmzp = 0.25_wp * &
[1]	240	( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
	241
[1320]	242	tend(k,j,i) = tend(k,j,i) &
[2037]	243	& + ( kmzp * ( ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1) &
	244	& + ( w(k,j,i) - w(k,j-1,i) ) * ddy &
	245	& ) * rho_air_zw(k) &
	246	& - ( -vsws(j,i) ) &
	247	& ) * ddzw(k) * drho_air(k)
[1]	248	ENDIF
	249
[102]	250	!
	251	!-- Vertical diffusion at the first gridpoint below the top boundary,
	252	!-- if the momentum flux at the top is prescribed by the user
[103]	253	IF ( use_top_fluxes .AND. constant_top_momentumflux ) THEN
[102]	254	k = nzt
	255	!
	256	!-- Interpolate eddy diffusivities on staggered gridpoints
[1340]	257	kmzm = 0.25_wp * &
[102]	258	( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
	259
[1320]	260	tend(k,j,i) = tend(k,j,i) &
[2037]	261	& + ( ( -vswst(j,i) ) &
	262	& - kmzm * ( ( v(k,j,i) - v(k-1,j,i) ) * ddzu(k) &
	263	& + ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy &
	264	& ) * rho_air_zw(k-1) &
	265	& ) * ddzw(k) * drho_air(k)
[102]	266	ENDIF
	267
[1]	268	ENDDO
	269	ENDDO
	270
	271	END SUBROUTINE diffusion_v
	272
	273
	274	!------------------------------------------------------------------------------!
[1682]	275	! Description:
	276	! ------------
	277	!> Call for all grid points - accelerator version
[1015]	278	!------------------------------------------------------------------------------!
	279	SUBROUTINE diffusion_v_acc
	280
[1320]	281	USE arrays_3d, &
[2037]	282	ONLY: ddzu, ddzw, km, tend, u, v, vsws, vswst, w, &
	283	drho_air, rho_air_zw
[1320]	284
	285	USE control_parameters, &
	286	ONLY: constant_top_momentumflux, topography, use_surface_fluxes, &
	287	use_top_fluxes
	288
	289	USE grid_variables, &
	290	ONLY: ddx, ddy, ddy2, fxm, fxp, wall_v
	291
	292	USE indices, &
	293	ONLY: i_left, i_right, j_north, j_south, nxl, nxr, nyn, nys, nzb, &
	294	nzb_diff_v, nzb_v_inner, nzb_v_outer, nzt, nzt_diff
	295
	296	USE kinds
[1015]	297
	298	IMPLICIT NONE
	299
[1682]	300	INTEGER(iwp) :: i !<
	301	INTEGER(iwp) :: j !<
	302	INTEGER(iwp) :: k !<
	303	REAL(wp) :: kmxm !<
	304	REAL(wp) :: kmxp !<
	305	REAL(wp) :: kmzm !<
	306	REAL(wp) :: kmzp !<
[1015]	307
[1682]	308	REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: vsus !<
[1015]	309	!$acc declare create ( vsus )
	310
	311	!
	312	!-- First calculate horizontal momentum flux v'u' at vertical walls,
	313	!-- if neccessary
	314	IF ( topography /= 'flat' ) THEN
[1320]	315	CALL wall_fluxes_acc( vsus, 0.0_wp, 1.0_wp, 0.0_wp, 0.0_wp, &
	316	nzb_v_inner, nzb_v_outer, wall_v )
[1015]	317	ENDIF
	318
[1320]	319	!$acc kernels present ( u, v, w, km, tend, vsws, vswst ) &
	320	!$acc present ( ddzu, ddzw, fxm, fxp, wall_v ) &
[1015]	321	!$acc present ( nzb_v_inner, nzb_v_outer, nzb_diff_v )
[1128]	322	DO i = i_left, i_right
	323	DO j = j_south, j_north
[1015]	324	!
	325	!-- Compute horizontal diffusion
	326	DO k = 1, nzt
	327	IF ( k > nzb_v_outer(j,i) ) THEN
	328	!
	329	!-- Interpolate eddy diffusivities on staggered gridpoints
[1340]	330	kmxp = 0.25_wp * &
[1015]	331	( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
[1340]	332	kmxm = 0.25_wp * &
[1015]	333	( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
	334
	335	tend(k,j,i) = tend(k,j,i) &
	336	& + ( kmxp * ( v(k,j,i+1) - v(k,j,i) ) * ddx &
	337	& + kmxp * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy &
	338	& - kmxm * ( v(k,j,i) - v(k,j,i-1) ) * ddx &
	339	& - kmxm * ( u(k,j,i) - u(k,j-1,i) ) * ddy &
	340	& ) * ddx &
[1340]	341	& + 2.0_wp * ( &
[1015]	342	& km(k,j,i) * ( v(k,j+1,i) - v(k,j,i) ) &
	343	& - km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) ) &
[1340]	344	& ) * ddy2
[1015]	345	ENDIF
	346	ENDDO
	347
	348	!
	349	!-- Wall functions at the left and right walls, respectively
	350	DO k = 1, nzt
[1320]	351	IF( k > nzb_v_inner(j,i) .AND. k <= nzb_v_outer(j,i) .AND. &
[1340]	352	wall_v(j,i) /= 0.0_wp ) THEN
[1015]	353
[1340]	354	kmxp = 0.25_wp * &
[1015]	355	( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
[1340]	356	kmxm = 0.25_wp * &
[1015]	357	( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
	358
	359	tend(k,j,i) = tend(k,j,i) &
[1340]	360	+ 2.0_wp * ( &
[1015]	361	km(k,j,i) * ( v(k,j+1,i) - v(k,j,i) ) &
	362	- km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) ) &
[1340]	363	) * ddy2 &
[1015]	364	+ ( fxp(j,i) * ( &
	365	kmxp * ( v(k,j,i+1) - v(k,j,i) ) * ddx &
	366	+ kmxp * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy &
	367	) &
	368	- fxm(j,i) * ( &
	369	kmxm * ( v(k,j,i) - v(k,j,i-1) ) * ddx &
	370	+ kmxm * ( u(k,j,i) - u(k,j-1,i) ) * ddy &
	371	) &
	372	+ wall_v(j,i) * vsus(k,j,i) &
	373	) * ddx
	374	ENDIF
	375	ENDDO
	376
	377	!
	378	!-- Compute vertical diffusion. In case of simulating a Prandtl
	379	!-- layer, index k starts at nzb_v_inner+2.
	380	DO k = 1, nzt_diff
	381	IF ( k >= nzb_diff_v(j,i) ) THEN
	382	!
	383	!-- Interpolate eddy diffusivities on staggered gridpoints
[1340]	384	kmzp = 0.25_wp * &
[1015]	385	( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
[1340]	386	kmzm = 0.25_wp * &
[1015]	387	( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
	388
	389	tend(k,j,i) = tend(k,j,i) &
	390	& + ( kmzp * ( ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1)&
	391	& + ( w(k,j,i) - w(k,j-1,i) ) * ddy &
[2037]	392	& ) * rho_air_zw(k) &
[1015]	393	& - kmzm * ( ( v(k,j,i) - v(k-1,j,i) ) * ddzu(k)&
	394	& + ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy &
[2037]	395	& ) * rho_air_zw(k-1) &
	396	& ) * ddzw(k) * drho_air(k)
[1015]	397	ENDIF
	398	ENDDO
	399
	400	ENDDO
	401	ENDDO
	402
	403	!
	404	!-- Vertical diffusion at the first grid point above the surface,
	405	!-- if the momentum flux at the bottom is given by the Prandtl law
	406	!-- or if it is prescribed by the user.
	407	!-- Difference quotient of the momentum flux is not formed over
	408	!-- half of the grid spacing (2.0*ddzw(k)) any more, since the
[1320]	409	!-- comparison with other (LES) models showed that the values of
[1015]	410	!-- the momentum flux becomes too large in this case.
	411	!-- The term containing w(k-1,..) (see above equation) is removed here
	412	!-- because the vertical velocity is assumed to be zero at the surface.
	413	IF ( use_surface_fluxes ) THEN
	414
[1128]	415	DO i = i_left, i_right
	416	DO j = j_south, j_north
[1015]	417
	418	k = nzb_v_inner(j,i)+1
	419	!
	420	!-- Interpolate eddy diffusivities on staggered gridpoints
[1340]	421	kmzp = 0.25_wp * &
[1015]	422	( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
	423
[1320]	424	tend(k,j,i) = tend(k,j,i) &
[2037]	425	& + ( kmzp * ( ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1) &
	426	& + ( w(k,j,i) - w(k,j-1,i) ) * ddy &
	427	& ) * rho_air_zw(k) &
	428	& - ( -vsws(j,i) ) &
	429	& ) * ddzw(k) * drho_air(k)
[1015]	430	ENDDO
	431	ENDDO
	432
	433	ENDIF
	434
	435	!
	436	!-- Vertical diffusion at the first gridpoint below the top boundary,
	437	!-- if the momentum flux at the top is prescribed by the user
	438	IF ( use_top_fluxes .AND. constant_top_momentumflux ) THEN
	439
	440	k = nzt
	441
[1128]	442	DO i = i_left, i_right
	443	DO j = j_south, j_north
[1015]	444
	445	!
	446	!-- Interpolate eddy diffusivities on staggered gridpoints
[1340]	447	kmzm = 0.25_wp * &
[1015]	448	( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
	449
[1320]	450	tend(k,j,i) = tend(k,j,i) &
[2037]	451	& + ( ( -vswst(j,i) ) &
	452	& - kmzm * ( ( v(k,j,i) - v(k-1,j,i) ) * ddzu(k) &
	453	& + ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy &
	454	& ) * rho_air_zw(k-1) &
	455	& ) * ddzw(k) * drho_air(k)
[1015]	456	ENDDO
	457	ENDDO
	458
	459	ENDIF
	460	!$acc end kernels
	461
	462	END SUBROUTINE diffusion_v_acc
	463
	464
	465	!------------------------------------------------------------------------------!
[1682]	466	! Description:
	467	! ------------
	468	!> Call for grid point i,j
[1]	469	!------------------------------------------------------------------------------!
[1001]	470	SUBROUTINE diffusion_v_ij( i, j )
[1]	471
[1320]	472	USE arrays_3d, &
[2037]	473	ONLY: ddzu, ddzw, km, tend, u, v, vsws, vswst, w, &
	474	drho_air, rho_air_zw
[1320]	475
	476	USE control_parameters, &
	477	ONLY: constant_top_momentumflux, use_surface_fluxes, use_top_fluxes
	478
	479	USE grid_variables, &
	480	ONLY: ddx, ddy, ddy2, fxm, fxp, wall_v
	481
	482	USE indices, &
	483	ONLY: nzb, nzb_diff_v, nzb_v_inner, nzb_v_outer, nzt, nzt_diff
	484
	485	USE kinds
[1]	486
	487	IMPLICIT NONE
	488
[1682]	489	INTEGER(iwp) :: i !<
	490	INTEGER(iwp) :: j !<
	491	INTEGER(iwp) :: k !<
	492	REAL(wp) :: kmxm !<
	493	REAL(wp) :: kmxp !<
	494	REAL(wp) :: kmzm !<
	495	REAL(wp) :: kmzp !<
[1]	496
[1682]	497	REAL(wp), DIMENSION(nzb:nzt+1) :: vsus !<
[1001]	498
[1]	499	!
	500	!-- Compute horizontal diffusion
	501	DO k = nzb_v_outer(j,i)+1, nzt
	502	!
	503	!-- Interpolate eddy diffusivities on staggered gridpoints
[1340]	504	kmxp = 0.25_wp * ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
	505	kmxm = 0.25_wp * ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
[1]	506
[1320]	507	tend(k,j,i) = tend(k,j,i) &
	508	& + ( kmxp * ( v(k,j,i+1) - v(k,j,i) ) * ddx &
	509	& + kmxp * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy &
	510	& - kmxm * ( v(k,j,i) - v(k,j,i-1) ) * ddx &
	511	& - kmxm * ( u(k,j,i) - u(k,j-1,i) ) * ddy &
	512	& ) * ddx &
[1340]	513	& + 2.0_wp * ( &
[1320]	514	& km(k,j,i) * ( v(k,j+1,i) - v(k,j,i) ) &
	515	& - km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) ) &
[1340]	516	& ) * ddy2
[1]	517	ENDDO
	518
	519	!
	520	!-- Wall functions at the left and right walls, respectively
[1340]	521	IF ( wall_v(j,i) /= 0.0_wp ) THEN
[51]	522
	523	!
	524	!-- Calculate the horizontal momentum flux v'u'
[1320]	525	CALL wall_fluxes( i, j, nzb_v_inner(j,i)+1, nzb_v_outer(j,i), &
	526	vsus, 0.0_wp, 1.0_wp, 0.0_wp, 0.0_wp )
[51]	527
[1]	528	DO k = nzb_v_inner(j,i)+1, nzb_v_outer(j,i)
[1340]	529	kmxp = 0.25_wp * &
[978]	530	( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
[1340]	531	kmxm = 0.25_wp * &
[978]	532	( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
[1]	533
	534	tend(k,j,i) = tend(k,j,i) &
[1340]	535	+ 2.0_wp * ( &
[1]	536	km(k,j,i) * ( v(k,j+1,i) - v(k,j,i) ) &
	537	- km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) ) &
[1340]	538	) * ddy2 &
[1]	539	+ ( fxp(j,i) * ( &
[978]	540	kmxp * ( v(k,j,i+1) - v(k,j,i) ) * ddx &
	541	+ kmxp * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy &
[1]	542	) &
	543	- fxm(j,i) * ( &
[978]	544	kmxm * ( v(k,j,i) - v(k,j,i-1) ) * ddx &
	545	+ kmxm * ( u(k,j,i) - u(k,j-1,i) ) * ddy &
[1]	546	) &
[51]	547	+ wall_v(j,i) * vsus(k) &
[1]	548	) * ddx
	549	ENDDO
	550	ENDIF
	551
	552	!
	553	!-- Compute vertical diffusion. In case of simulating a Prandtl layer,
	554	!-- index k starts at nzb_v_inner+2.
[102]	555	DO k = nzb_diff_v(j,i), nzt_diff
[1]	556	!
	557	!-- Interpolate eddy diffusivities on staggered gridpoints
[1340]	558	kmzp = 0.25_wp * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
	559	kmzm = 0.25_wp * ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
[1]	560
[1320]	561	tend(k,j,i) = tend(k,j,i) &
	562	& + ( kmzp * ( ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1) &
	563	& + ( w(k,j,i) - w(k,j-1,i) ) * ddy &
[2037]	564	& ) * rho_air_zw(k) &
[1320]	565	& - kmzm * ( ( v(k,j,i) - v(k-1,j,i) ) * ddzu(k) &
	566	& + ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy &
[2037]	567	& ) * rho_air_zw(k-1) &
	568	& ) * ddzw(k) * drho_air(k)
[1]	569	ENDDO
	570
	571	!
	572	!-- Vertical diffusion at the first grid point above the surface, if the
	573	!-- momentum flux at the bottom is given by the Prandtl law or if it is
	574	!-- prescribed by the user.
	575	!-- Difference quotient of the momentum flux is not formed over half of
	576	!-- the grid spacing (2.0*ddzw(k)) any more, since the comparison with
[1320]	577	!-- other (LES) models showed that the values of the momentum flux becomes
[1]	578	!-- too large in this case.
	579	!-- The term containing w(k-1,..) (see above equation) is removed here
	580	!-- because the vertical velocity is assumed to be zero at the surface.
	581	IF ( use_surface_fluxes ) THEN
	582	k = nzb_v_inner(j,i)+1
	583	!
	584	!-- Interpolate eddy diffusivities on staggered gridpoints
[1340]	585	kmzp = 0.25_wp * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
[1]	586
[1320]	587	tend(k,j,i) = tend(k,j,i) &
[2037]	588	& + ( kmzp * ( ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1) &
	589	& + ( w(k,j,i) - w(k,j-1,i) ) * ddy &
	590	& ) * rho_air_zw(k) &
	591	& - ( -vsws(j,i) ) &
	592	& ) * ddzw(k) * drho_air(k)
[1]	593	ENDIF
	594
[102]	595	!
	596	!-- Vertical diffusion at the first gridpoint below the top boundary,
	597	!-- if the momentum flux at the top is prescribed by the user
[103]	598	IF ( use_top_fluxes .AND. constant_top_momentumflux ) THEN
[102]	599	k = nzt
	600	!
	601	!-- Interpolate eddy diffusivities on staggered gridpoints
[2037]	602	kmzm = 0.25_wp * ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
[102]	603
[1320]	604	tend(k,j,i) = tend(k,j,i) &
[2037]	605	& + ( ( -vswst(j,i) ) &
	606	& - kmzm * ( ( v(k,j,i) - v(k-1,j,i) ) * ddzu(k) &
	607	& + ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy &
	608	& ) * rho_air_zw(k-1) &
	609	& ) * ddzw(k) * drho_air(k)
[102]	610	ENDIF
	611
[1]	612	END SUBROUTINE diffusion_v_ij
	613
[1321]	614	END MODULE diffusion_v_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |