Home

Context Navigation

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

Last change on this file since 1324 was 1321, checked in by raasch, 11 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 27.1 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |