Home

Context Navigation

diffusion_u_mod.f90 @ 1851

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

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

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

Context Navigation

source: palm/trunk/SOURCE/diffusion_u_mod.f90 @ 1851

Download in other formats: