Home

Context Navigation

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

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |