Home

Context Navigation

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

Last change on this file since 4775 was 4674, checked in by pavelkrc, 4 years ago
Update ACC directives for downward facing USM and LSM surfaces
Property svn:keywords set to `Id`
File size: 25.5 KB

Line
1	!> @file diffusion_u.f90
2	!--------------------------------------------------------------------------------------------------!
3	! This file is part of the PALM model system.
4	!
5	! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General
6	! Public License as published by the Free Software Foundation, either version 3 of the License, or
7	! (at your option) any later version.
8	!
9	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
10	! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
11	! Public License for more details.
12	!
13	! You should have received a copy of the GNU General Public License along with PALM. If not, see
14	! <http://www.gnu.org/licenses/>.
15	!
16	! Copyright 1997-2020 Leibniz Universitaet Hannover
17	!--------------------------------------------------------------------------------------------------!
18	!
19	! Current revisions:
20	! -----------------
21	!
22	!
23	! Former revisions:
24	! -----------------
25	! $Id: diffusion_u.f90 4674 2020-09-10 10:36:55Z raasch $
26	! Update ACC directives for downward facing USM and LSM surfaces
27	!
28	! 4671 2020-09-09 20:27:58Z pavelkrc
29	! Implementation of downward facing USM and LSM surfaces
30	!
31	! 4583 2020-06-29 12:36:47Z raasch
32	! file re-formatted to follow the PALM coding standard
33	!
34	! 4360 2020-01-07 11:25:50Z suehring
35	! Introduction of wall_flags_total_0, which currently sets bits based on static topography
36	! information used in wall_flags_static_0
37	!
38	! 4329 2019-12-10 15:46:36Z motisi
39	! Renamed wall_flags_0 to wall_flags_static_0
40	!
41	! 4182 2019-08-22 15:20:23Z scharf
42	! Corrected "Former revisions" section
43	!
44	! 3655 2019-01-07 16:51:22Z knoop
45	! OpenACC port for SPEC
46	!
47	! Revision 1.1 1997/09/12 06:23:51 raasch
48	! Initial revision
49	!
50	!
51	! Description:
52	! ------------
53	!> Diffusion term of the u-component
54	!> @todo additional damping (needed for non-cyclic bc) causes bad vectorization and slows down the
55	! speed on NEC about 5-10%
56	!--------------------------------------------------------------------------------------------------!
57	MODULE diffusion_u_mod
58
59
60	PRIVATE
61	PUBLIC diffusion_u
62
63	INTERFACE diffusion_u
64	MODULE PROCEDURE diffusion_u
65	MODULE PROCEDURE diffusion_u_ij
66	END INTERFACE diffusion_u
67
68	CONTAINS
69
70
71	!--------------------------------------------------------------------------------------------------!
72	! Description:
73	! ------------
74	!> Call for all grid points
75	!--------------------------------------------------------------------------------------------------!
76	SUBROUTINE diffusion_u
77
78	USE arrays_3d, &
79	ONLY: ddzu, ddzw, drho_air, km, rho_air_zw, tend, u, v, w
80
81	USE control_parameters, &
82	ONLY: constant_top_momentumflux, use_surface_fluxes, use_top_fluxes
83
84	USE grid_variables, &
85	ONLY: ddx, ddx2, ddy
86
87	USE indices, &
88	ONLY: nxlu, nxr, nyn, nys, nzb, nzt, wall_flags_total_0
89
90	USE kinds
91
92	USE surface_mod, &
93	ONLY : surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, surf_usm_v
94
95	IMPLICIT NONE
96
97	INTEGER(iwp) :: i !< running index x direction
98	INTEGER(iwp) :: j !< running index y direction
99	INTEGER(iwp) :: k !< running index z direction
100	INTEGER(iwp) :: l !< running index of surface type, south- or north-facing wall
101	INTEGER(iwp) :: m !< running index surface elements
102	INTEGER(iwp) :: surf_e !< end index of surface elements at (j,i)-gridpoint
103	INTEGER(iwp) :: surf_s !< start index of surface elements at (j,i)-gridpoint
104
105	REAL(wp) :: flag !< flag to mask topography grid points
106	REAL(wp) :: kmym !< diffusion coefficient on southward side of the u-gridbox - interpolated onto xu-yv grid
107	REAL(wp) :: kmyp !< diffusion coefficient on northward side of the u-gridbox - interpolated onto xu-yv grid
108	REAL(wp) :: kmzm !< diffusion coefficient on bottom of the gridbox - interpolated onto xu-zw grid
109	REAL(wp) :: kmzp !< diffusion coefficient on top of the gridbox - interpolated onto xu-zw grid
110	REAL(wp) :: mask_bottom !< flag to mask vertical upward-facing surface
111	REAL(wp) :: mask_north !< flag to mask vertical surface north of the grid point
112	REAL(wp) :: mask_south !< flag to mask vertical surface south of the grid point
113	REAL(wp) :: mask_top !< flag to mask vertical downward-facing surface
114
115
116
117	!$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j, k, l, m) &
118	!$ACC PRIVATE(surf_e, surf_s, flag, kmym, kmyp, kmzm, kmzp) &
119	!$ACC PRIVATE(mask_bottom, mask_north, mask_south, mask_top) &
120	!$ACC PRESENT(wall_flags_total_0, km) &
121	!$ACC PRESENT(u, v, w) &
122	!$ACC PRESENT(ddzu, ddzw, drho_air, rho_air_zw) &
123	!$ACC PRESENT(surf_def_h(0:2), surf_def_v(0:1)) &
124	!$ACC PRESENT(surf_lsm_h(0:1), surf_lsm_v(0:1)) &
125	!$ACC PRESENT(surf_usm_h(0:1), surf_usm_v(0:1)) &
126	!$ACC PRESENT(tend)
127	DO i = nxlu, nxr
128	DO j = nys, nyn
129	!
130	!-- Compute horizontal diffusion
131	DO k = nzb+1, nzt
132	!
133	!-- Predetermine flag to mask topography and wall-bounded grid points.
134	!-- It is sufficient to masked only north- and south-facing surfaces, which need special
135	!-- treatment for the u-component.
136	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 1 ) )
137	mask_south = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j-1,i), 1 ) )
138	mask_north = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j+1,i), 1 ) )
139	!
140	!-- Interpolate eddy diffusivities on staggered gridpoints
141	kmyp = 0.25_wp * ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
142	kmym = 0.25_wp * ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
143
144	tend(k,j,i) = tend(k,j,i) &
145	+ 2.0_wp * ( &
146	km(k,j,i) * ( u(k,j,i+1) - u(k,j,i) ) &
147	- km(k,j,i-1) * ( u(k,j,i) - u(k,j,i-1) ) &
148	) * ddx2 * flag &
149	+ ( mask_north * ( &
150	kmyp * ( u(k,j+1,i) - u(k,j,i) ) * ddy &
151	+ kmyp * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx &
152	) &
153	- mask_south * ( &
154	kmym * ( u(k,j,i) - u(k,j-1,i) ) * ddy &
155	+ kmym * ( v(k,j,i) - v(k,j,i-1) ) * ddx &
156	) &
157	) * ddy * flag
158	ENDDO
159	!
160	!-- Add horizontal momentum flux u'v' at north- (l=0) and south-facing (l=1) surfaces.
161	!-- Note, in the the flat case, loops won't be entered as start_index > end_index.
162	!-- Furtermore, note, no vertical natural surfaces so far.
163	!-- Default-type surfaces
164	DO l = 0, 1
165	surf_s = surf_def_v(l)%start_index(j,i)
166	surf_e = surf_def_v(l)%end_index(j,i)
167	DO m = surf_s, surf_e
168	k = surf_def_v(l)%k(m)
169	tend(k,j,i) = tend(k,j,i) + surf_def_v(l)%mom_flux_uv(m) * ddy
170	ENDDO
171	ENDDO
172	!
173	!-- Natural-type surfaces
174	DO l = 0, 1
175	surf_s = surf_lsm_v(l)%start_index(j,i)
176	surf_e = surf_lsm_v(l)%end_index(j,i)
177	DO m = surf_s, surf_e
178	k = surf_lsm_v(l)%k(m)
179	tend(k,j,i) = tend(k,j,i) + surf_lsm_v(l)%mom_flux_uv(m) * ddy
180	ENDDO
181	ENDDO
182	!
183	!-- Urban-type surfaces
184	DO l = 0, 1
185	surf_s = surf_usm_v(l)%start_index(j,i)
186	surf_e = surf_usm_v(l)%end_index(j,i)
187	DO m = surf_s, surf_e
188	k = surf_usm_v(l)%k(m)
189	tend(k,j,i) = tend(k,j,i) + surf_usm_v(l)%mom_flux_uv(m) * ddy
190	ENDDO
191	ENDDO
192
193	!
194	!-- Compute vertical diffusion. In case of simulating a surface layer, respective grid
195	!-- diffusive fluxes are masked (flag 8) within this loop, and added further below, else,
196	!-- simple gradient approach is applied. Model top is also mask if top-momentum flux is
197	!-- given.
198	DO k = nzb+1, nzt
199	!
200	!-- Determine flags to mask topography below and above. Flag 1 is used to mask
201	!-- topography in general, and flag 8 implies information about use_surface_fluxes.
202	!-- Flag 9 is used to control momentum flux at model top.
203	mask_bottom = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k-1,j,i), 8 ) )
204	mask_top = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k+1,j,i), 8 ) ) * &
205	MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k+1,j,i), 9 ) )
206	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 1 ) )
207	!
208	!-- Interpolate eddy diffusivities on staggered gridpoints
209	kmzp = 0.25_wp * ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
210	kmzm = 0.25_wp * ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
211
212	tend(k,j,i) = tend(k,j,i) &
213	+ ( kmzp * ( ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1) &
214	+ ( w(k,j,i) - w(k,j,i-1) ) * ddx &
215	) * rho_air_zw(k) * mask_top &
216	- kmzm * ( ( u(k,j,i) - u(k-1,j,i) ) * ddzu(k) &
217	+ ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx &
218	) * rho_air_zw(k-1) * mask_bottom &
219	) * ddzw(k) * drho_air(k) * flag
220	ENDDO
221
222	!
223	!-- Vertical diffusion at the first grid point above the surface, if the momentum flux at
224	!-- the bottom is given by the Prandtl law or if it is prescribed by the user.
225	!-- Difference quotient of the momentum flux is not formed over half of the grid spacing
226	!-- (2.0*ddzw(k)) any more, since the comparison with other (LES) models showed that the
227	!-- values of the momentum flux becomes too large in this case.
228	!-- The term containing w(k-1,..) (see above equation) is removed here because the vertical
229	!-- velocity is assumed to be zero at the surface.
230	IF ( use_surface_fluxes ) THEN
231	!
232	!-- Default-type surfaces, upward-facing
233	surf_s = surf_def_h(0)%start_index(j,i)
234	surf_e = surf_def_h(0)%end_index(j,i)
235	DO m = surf_s, surf_e
236
237	k = surf_def_h(0)%k(m)
238
239	tend(k,j,i) = tend(k,j,i) &
240	+ ( - ( - surf_def_h(0)%usws(m) ) ) * ddzw(k) * drho_air(k)
241	ENDDO
242	!
243	!-- Default-type surfaces, dowward-facing
244	surf_s = surf_def_h(1)%start_index(j,i)
245	surf_e = surf_def_h(1)%end_index(j,i)
246	DO m = surf_s, surf_e
247
248	k = surf_def_h(1)%k(m)
249
250	tend(k,j,i) = tend(k,j,i) &
251	+ ( - surf_def_h(1)%usws(m) ) * ddzw(k) * drho_air(k)
252	ENDDO
253	!
254	!-- Natural-type surfaces, upward-facing
255	surf_s = surf_lsm_h(0)%start_index(j,i)
256	surf_e = surf_lsm_h(0)%end_index(j,i)
257	DO m = surf_s, surf_e
258
259	k = surf_lsm_h(0)%k(m)
260
261	tend(k,j,i) = tend(k,j,i) &
262	+ ( - ( - surf_lsm_h(0)%usws(m) ) ) * ddzw(k) * drho_air(k)
263	ENDDO
264	!
265	!-- Natural-type surfaces, downward-facing
266	surf_s = surf_lsm_h(1)%start_index(j,i)
267	surf_e = surf_lsm_h(1)%end_index(j,i)
268	DO m = surf_s, surf_e
269
270	k = surf_lsm_h(1)%k(m)
271
272	tend(k,j,i) = tend(k,j,i) &
273	+ ( - surf_lsm_h(1)%usws(m) ) * ddzw(k) * drho_air(k)
274	ENDDO
275	!
276	!-- Urban-type surfaces, upward-facing
277	surf_s = surf_usm_h(0)%start_index(j,i)
278	surf_e = surf_usm_h(0)%end_index(j,i)
279	DO m = surf_s, surf_e
280
281	k = surf_usm_h(0)%k(m)
282
283	tend(k,j,i) = tend(k,j,i) &
284	+ ( - ( - surf_usm_h(0)%usws(m) ) ) * ddzw(k) * drho_air(k)
285	ENDDO
286	!
287	!-- Urban-type surfaces, downward-facing
288	surf_s = surf_usm_h(1)%start_index(j,i)
289	surf_e = surf_usm_h(1)%end_index(j,i)
290	DO m = surf_s, surf_e
291
292	k = surf_usm_h(1)%k(m)
293
294	tend(k,j,i) = tend(k,j,i) &
295	+ ( - surf_usm_h(1)%usws(m) ) * ddzw(k) * drho_air(k)
296	ENDDO
297
298	ENDIF
299	!
300	!-- Add momentum flux at model top
301	IF ( use_top_fluxes .AND. constant_top_momentumflux ) THEN
302	surf_s = surf_def_h(2)%start_index(j,i)
303	surf_e = surf_def_h(2)%end_index(j,i)
304	DO m = surf_s, surf_e
305
306	k = surf_def_h(2)%k(m)
307
308	tend(k,j,i) = tend(k,j,i) &
309	+ ( - surf_def_h(2)%usws(m) ) * ddzw(k) * drho_air(k)
310	ENDDO
311	ENDIF
312
313	ENDDO
314	ENDDO
315
316	END SUBROUTINE diffusion_u
317
318
319	!--------------------------------------------------------------------------------------------------!
320	! Description:
321	! ------------
322	!> Call for grid point i,j
323	!--------------------------------------------------------------------------------------------------!
324	SUBROUTINE diffusion_u_ij( i, j )
325
326	USE arrays_3d, &
327	ONLY: ddzu, ddzw, drho_air, km, tend, u, v, w, rho_air_zw
328
329	USE control_parameters, &
330	ONLY: constant_top_momentumflux, use_surface_fluxes, use_top_fluxes
331
332	USE grid_variables, &
333	ONLY: ddx, ddx2, ddy
334
335	USE indices, &
336	ONLY: nzb, nzt, wall_flags_total_0
337
338	USE kinds
339
340	USE surface_mod, &
341	ONLY : surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, surf_usm_v
342
343	IMPLICIT NONE
344
345	INTEGER(iwp) :: i !< running index x direction
346	INTEGER(iwp) :: j !< running index y direction
347	INTEGER(iwp) :: k !< running index z direction
348	INTEGER(iwp) :: l !< running index of surface type, south- or north-facing wall
349	INTEGER(iwp) :: m !< running index surface elements
350	INTEGER(iwp) :: surf_e !< End index of surface elements at (j,i)-gridpoint
351	INTEGER(iwp) :: surf_s !< Start index of surface elements at (j,i)-gridpoint
352
353	REAL(wp) :: flag !< flag to mask topography grid points
354	REAL(wp) :: kmym !< diffusion coefficient on southward side of the u-gridbox - interpolated onto xu-yv grid
355	REAL(wp) :: kmyp !<diffusion coefficient on northward side of the u-gridbox - interpolated onto xu-yv grid
356	REAL(wp) :: kmzm !< diffusion coefficient on bottom of the gridbox - interpolated onto xu-zw grid
357	REAL(wp) :: kmzp !< diffusion coefficient on top of the gridbox - interpolated onto xu-zw grid
358	REAL(wp) :: mask_bottom !< flag to mask vertical upward-facing surface
359	REAL(wp) :: mask_north !< flag to mask vertical surface north of the grid point
360	REAL(wp) :: mask_south !< flag to mask vertical surface south of the grid point
361	REAL(wp) :: mask_top !< flag to mask vertical downward-facing surface
362	!
363	!-- Compute horizontal diffusion
364	DO k = nzb+1, nzt
365	!
366	!-- Predetermine flag to mask topography and wall-bounded grid points.
367	!-- It is sufficient to masked only north- and south-facing surfaces, which need special
368	!-- treatment for the u-component.
369	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 1 ) )
370	mask_south = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j-1,i), 1 ) )
371	mask_north = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j+1,i), 1 ) )
372	!
373	!-- Interpolate eddy diffusivities on staggered gridpoints
374	kmyp = 0.25_wp * ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
375	kmym = 0.25_wp * ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
376
377	tend(k,j,i) = tend(k,j,i) &
378	+ 2.0_wp * ( &
379	km(k,j,i) * ( u(k,j,i+1) - u(k,j,i) ) &
380	- km(k,j,i-1) * ( u(k,j,i) - u(k,j,i-1) ) &
381	) * ddx2 * flag &
382	+ ( &
383	mask_north * kmyp * ( ( u(k,j+1,i) - u(k,j,i) ) * ddy &
384	+ ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx &
385	) &
386	- mask_south * kmym * ( ( u(k,j,i) - u(k,j-1,i) ) * ddy &
387	+ ( v(k,j,i) - v(k,j,i-1) ) * ddx &
388	) &
389	) * ddy * flag
390	ENDDO
391
392	!
393	!-- Add horizontal momentum flux u'v' at north- (l=0) and south-facing (l=1) surfaces. Note, in
394	!-- the the flat case, loops won't be entered as start_index > end_index. Furtermore, note, no
395	!-- vertical natural surfaces so far.
396	!-- Default-type surfaces
397	DO l = 0, 1
398	surf_s = surf_def_v(l)%start_index(j,i)
399	surf_e = surf_def_v(l)%end_index(j,i)
400	DO m = surf_s, surf_e
401	k = surf_def_v(l)%k(m)
402	tend(k,j,i) = tend(k,j,i) + surf_def_v(l)%mom_flux_uv(m) * ddy
403	ENDDO
404	ENDDO
405	!
406	!-- Natural-type surfaces
407	DO l = 0, 1
408	surf_s = surf_lsm_v(l)%start_index(j,i)
409	surf_e = surf_lsm_v(l)%end_index(j,i)
410	DO m = surf_s, surf_e
411	k = surf_lsm_v(l)%k(m)
412	tend(k,j,i) = tend(k,j,i) + surf_lsm_v(l)%mom_flux_uv(m) * ddy
413	ENDDO
414	ENDDO
415	!
416	!-- Urban-type surfaces
417	DO l = 0, 1
418	surf_s = surf_usm_v(l)%start_index(j,i)
419	surf_e = surf_usm_v(l)%end_index(j,i)
420	DO m = surf_s, surf_e
421	k = surf_usm_v(l)%k(m)
422	tend(k,j,i) = tend(k,j,i) + surf_usm_v(l)%mom_flux_uv(m) * ddy
423	ENDDO
424	ENDDO
425	!
426	!-- Compute vertical diffusion. In case of simulating a surface layer, respective grid diffusive
427	!-- fluxes are masked (flag 8) within this loop, and added further below, else, simple gradient
428	!-- approach is applied. Model top is also mask if top-momentum flux is given.
429	DO k = nzb+1, nzt
430	!
431	!-- Determine flags to mask topography below and above. Flag 1 is used to mask topography in
432	!-- general, and flag 8 implies information about use_surface_fluxes. Flag 9 is used to
433	!-- control momentum flux at model top.
434	mask_bottom = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k-1,j,i), 8 ) )
435	mask_top = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k+1,j,i), 8 ) ) * &
436	MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k+1,j,i), 9 ) )
437	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 1 ) )
438	!
439	!-- Interpolate eddy diffusivities on staggered gridpoints
440	kmzp = 0.25_wp * ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
441	kmzm = 0.25_wp * ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
442
443	tend(k,j,i) = tend(k,j,i) &
444	+ ( kmzp * ( ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1) &
445	+ ( w(k,j,i) - w(k,j,i-1) ) * ddx &
446	) * rho_air_zw(k) * mask_top &
447	- kmzm * ( ( u(k,j,i) - u(k-1,j,i) ) * ddzu(k) &
448	+ ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx &
449	) * rho_air_zw(k-1) * mask_bottom &
450	) * ddzw(k) * drho_air(k) * flag
451	ENDDO
452
453	!
454	!-- Vertical diffusion at the first surface grid points, if the momentum flux at the bottom is
455	!-- given by the Prandtl law or if it is prescribed by the user.
456	!-- Difference quotient of the momentum flux is not formed over half of the grid spacing
457	!-- (2.0*ddzw(k)) any more, since the comparison with other (LES) models showed that the values
458	!-- of the momentum flux becomes too large in this case.
459	IF ( use_surface_fluxes ) THEN
460	!
461	!-- Default-type surfaces, upward-facing
462	surf_s = surf_def_h(0)%start_index(j,i)
463	surf_e = surf_def_h(0)%end_index(j,i)
464	DO m = surf_s, surf_e
465
466	k = surf_def_h(0)%k(m)
467
468	tend(k,j,i) = tend(k,j,i) + ( - ( - surf_def_h(0)%usws(m) ) ) * ddzw(k) * drho_air(k)
469	ENDDO
470	!
471	!-- Default-type surfaces, dowward-facing (except for model-top fluxes)
472	surf_s = surf_def_h(1)%start_index(j,i)
473	surf_e = surf_def_h(1)%end_index(j,i)
474	DO m = surf_s, surf_e
475
476	k = surf_def_h(1)%k(m)
477
478	tend(k,j,i) = tend(k,j,i) + ( - surf_def_h(1)%usws(m) ) * ddzw(k) * drho_air(k)
479	ENDDO
480	!
481	!-- Natural-type surfaces, upward-facing
482	surf_s = surf_lsm_h(0)%start_index(j,i)
483	surf_e = surf_lsm_h(0)%end_index(j,i)
484	DO m = surf_s, surf_e
485
486	k = surf_lsm_h(0)%k(m)
487
488	tend(k,j,i) = tend(k,j,i) + ( - ( - surf_lsm_h(0)%usws(m) ) ) * ddzw(k) * drho_air(k)
489	ENDDO
490	!
491	!-- Natural-type surfaces, downward-facing
492	surf_s = surf_lsm_h(1)%start_index(j,i)
493	surf_e = surf_lsm_h(1)%end_index(j,i)
494	DO m = surf_s, surf_e
495
496	k = surf_lsm_h(1)%k(m)
497
498	tend(k,j,i) = tend(k,j,i) + ( - surf_lsm_h(1)%usws(m) ) * ddzw(k) * drho_air(k)
499	ENDDO
500	!
501	!-- Urban-type surfaces, upward-facing
502	surf_s = surf_usm_h(0)%start_index(j,i)
503	surf_e = surf_usm_h(0)%end_index(j,i)
504	DO m = surf_s, surf_e
505
506	k = surf_usm_h(0)%k(m)
507
508	tend(k,j,i) = tend(k,j,i) + ( - ( - surf_usm_h(0)%usws(m) ) ) * ddzw(k) * drho_air(k)
509	ENDDO
510	!
511	!-- Urban-type surfaces, downward-facing
512	surf_s = surf_usm_h(1)%start_index(j,i)
513	surf_e = surf_usm_h(1)%end_index(j,i)
514	DO m = surf_s, surf_e
515
516	k = surf_usm_h(1)%k(m)
517
518	tend(k,j,i) = tend(k,j,i) + ( - surf_usm_h(1)%usws(m) ) * ddzw(k) * drho_air(k)
519	ENDDO
520
521	ENDIF
522	!
523	!-- Add momentum flux at model top
524	IF ( use_top_fluxes .AND. constant_top_momentumflux ) THEN
525	surf_s = surf_def_h(2)%start_index(j,i)
526	surf_e = surf_def_h(2)%end_index(j,i)
527	DO m = surf_s, surf_e
528
529	k = surf_def_h(2)%k(m)
530
531	tend(k,j,i) = tend(k,j,i) + ( - surf_def_h(2)%usws(m) ) * ddzw(k) * drho_air(k)
532	ENDDO
533	ENDIF
534
535
536	END SUBROUTINE diffusion_u_ij
537
538	END MODULE diffusion_u_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |