Home

Context Navigation

source: palm/trunk/SOURCE/diffusion_s.f90 @ 4331

Last change on this file since 4331 was 4329, checked in by motisi, 5 years ago
Renamed wall_flags_0 to wall_flags_static_0
Property svn:keywords set to `Id`
File size: 32.9 KB

Line
1	!> @file diffusion_s.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
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.
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	!
17	! Copyright 1997-2019 Leibniz Universitaet Hannover
18	!------------------------------------------------------------------------------!
19	!
20	! Current revisions:
21	! ------------------
22	!
23	!
24	! Former revisions:
25	! -----------------
26	! $Id: diffusion_s.f90 4329 2019-12-10 15:46:36Z suehring $
27	! Renamed wall_flags_0 to wall_flags_static_0
28	!
29	! 4182 2019-08-22 15:20:23Z scharf
30	! Corrected "Former revisions" section
31	!
32	! 3927 2019-04-23 13:24:29Z raasch
33	! pointer attribute removed from scalar 3d-array for performance reasons
34	!
35	! 3761 2019-02-25 15:31:42Z raasch
36	! unused variables removed
37	!
38	! 3655 2019-01-07 16:51:22Z knoop
39	! nopointer option removed
40	!
41	! Revision 1.1 2000/04/13 14:54:02 schroeter
42	! Initial revision
43	!
44	!
45	! Description:
46	! ------------
47	!> Diffusion term of scalar quantities (temperature and water content)
48	!------------------------------------------------------------------------------!
49	MODULE diffusion_s_mod
50
51
52	PRIVATE
53	PUBLIC diffusion_s
54
55	INTERFACE diffusion_s
56	MODULE PROCEDURE diffusion_s
57	MODULE PROCEDURE diffusion_s_ij
58	END INTERFACE diffusion_s
59
60	CONTAINS
61
62
63	!------------------------------------------------------------------------------!
64	! Description:
65	! ------------
66	!> Call for all grid points
67	!------------------------------------------------------------------------------!
68	SUBROUTINE diffusion_s( s, s_flux_def_h_up, s_flux_def_h_down, &
69	s_flux_t, &
70	s_flux_lsm_h_up, s_flux_usm_h_up, &
71	s_flux_def_v_north, s_flux_def_v_south, &
72	s_flux_def_v_east, s_flux_def_v_west, &
73	s_flux_lsm_v_north, s_flux_lsm_v_south, &
74	s_flux_lsm_v_east, s_flux_lsm_v_west, &
75	s_flux_usm_v_north, s_flux_usm_v_south, &
76	s_flux_usm_v_east, s_flux_usm_v_west )
77
78	USE arrays_3d, &
79	ONLY: ddzu, ddzw, kh, tend, drho_air, rho_air_zw
80
81	USE control_parameters, &
82	ONLY: use_surface_fluxes, use_top_fluxes
83
84	USE grid_variables, &
85	ONLY: ddx, ddx2, ddy, ddy2
86
87	USE indices, &
88	ONLY: nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt, wall_flags_static_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, &
94	surf_usm_v
95
96	IMPLICIT NONE
97
98	INTEGER(iwp) :: i !< running index x direction
99	INTEGER(iwp) :: j !< running index y direction
100	INTEGER(iwp) :: k !< running index z direction
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) :: mask_bottom !< flag to mask vertical upward-facing surface
107	REAL(wp) :: mask_east !< flag to mask vertical surface east of the grid point
108	REAL(wp) :: mask_north !< flag to mask vertical surface north of the grid point
109	REAL(wp) :: mask_south !< flag to mask vertical surface south of the grid point
110	REAL(wp) :: mask_west !< flag to mask vertical surface west of the grid point
111	REAL(wp) :: mask_top !< flag to mask vertical downward-facing surface
112
113	REAL(wp), DIMENSION(1:surf_def_v(0)%ns) :: s_flux_def_v_north !< flux at north-facing vertical default-type surfaces
114	REAL(wp), DIMENSION(1:surf_def_v(1)%ns) :: s_flux_def_v_south !< flux at south-facing vertical default-type surfaces
115	REAL(wp), DIMENSION(1:surf_def_v(2)%ns) :: s_flux_def_v_east !< flux at east-facing vertical default-type surfaces
116	REAL(wp), DIMENSION(1:surf_def_v(3)%ns) :: s_flux_def_v_west !< flux at west-facing vertical default-type surfaces
117	REAL(wp), DIMENSION(1:surf_def_h(0)%ns) :: s_flux_def_h_up !< flux at horizontal upward-facing default-type surfaces
118	REAL(wp), DIMENSION(1:surf_def_h(1)%ns) :: s_flux_def_h_down !< flux at horizontal donwward-facing default-type surfaces
119	REAL(wp), DIMENSION(1:surf_lsm_h%ns) :: s_flux_lsm_h_up !< flux at horizontal upward-facing natural-type surfaces
120	REAL(wp), DIMENSION(1:surf_lsm_v(0)%ns) :: s_flux_lsm_v_north !< flux at north-facing vertical natural-type surfaces
121	REAL(wp), DIMENSION(1:surf_lsm_v(1)%ns) :: s_flux_lsm_v_south !< flux at south-facing vertical natural-type surfaces
122	REAL(wp), DIMENSION(1:surf_lsm_v(2)%ns) :: s_flux_lsm_v_east !< flux at east-facing vertical natural-type surfaces
123	REAL(wp), DIMENSION(1:surf_lsm_v(3)%ns) :: s_flux_lsm_v_west !< flux at west-facing vertical natural-type surfaces
124	REAL(wp), DIMENSION(1:surf_usm_h%ns) :: s_flux_usm_h_up !< flux at horizontal upward-facing urban-type surfaces
125	REAL(wp), DIMENSION(1:surf_usm_v(0)%ns) :: s_flux_usm_v_north !< flux at north-facing vertical urban-type surfaces
126	REAL(wp), DIMENSION(1:surf_usm_v(1)%ns) :: s_flux_usm_v_south !< flux at south-facing vertical urban-type surfaces
127	REAL(wp), DIMENSION(1:surf_usm_v(2)%ns) :: s_flux_usm_v_east !< flux at east-facing vertical urban-type surfaces
128	REAL(wp), DIMENSION(1:surf_usm_v(3)%ns) :: s_flux_usm_v_west !< flux at west-facing vertical urban-type surfaces
129	REAL(wp), DIMENSION(1:surf_def_h(2)%ns) :: s_flux_t !< flux at model top
130
131	REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: s !< treated scalar
132
133
134	!$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j, k, m) &
135	!$ACC PRIVATE(surf_e, surf_s, flag, mask_top, mask_bottom) &
136	!$ACC PRIVATE(mask_north, mask_south, mask_west, mask_east) &
137	!$ACC PRESENT(wall_flags_static_0, kh) &
138	!$ACC PRESENT(s) &
139	!$ACC PRESENT(ddzu, ddzw, drho_air, rho_air_zw) &
140	!$ACC PRESENT(surf_def_h(0:2), surf_def_v(0:3)) &
141	!$ACC PRESENT(surf_lsm_h, surf_lsm_v(0:3)) &
142	!$ACC PRESENT(surf_usm_h, surf_usm_v(0:3)) &
143	!$ACC PRESENT(s_flux_def_h_up, s_flux_def_h_down) &
144	!$ACC PRESENT(s_flux_t) &
145	!$ACC PRESENT(s_flux_def_v_north, s_flux_def_v_south) &
146	!$ACC PRESENT(s_flux_def_v_east, s_flux_def_v_west) &
147	!$ACC PRESENT(s_flux_lsm_h_up) &
148	!$ACC PRESENT(s_flux_lsm_v_north, s_flux_lsm_v_south) &
149	!$ACC PRESENT(s_flux_lsm_v_east, s_flux_lsm_v_west) &
150	!$ACC PRESENT(s_flux_usm_h_up) &
151	!$ACC PRESENT(s_flux_usm_v_north, s_flux_usm_v_south) &
152	!$ACC PRESENT(s_flux_usm_v_east, s_flux_usm_v_west) &
153	!$ACC PRESENT(tend)
154	DO i = nxl, nxr
155	DO j = nys,nyn
156	!
157	!-- Compute horizontal diffusion
158	DO k = nzb+1, nzt
159	!
160	!-- Predetermine flag to mask topography and wall-bounded grid points
161	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_static_0(k,j,i), 0 ) )
162	!
163	!-- Predetermine flag to mask wall-bounded grid points, equivalent to
164	!-- former s_outer array
165	mask_west = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_static_0(k,j,i-1), 0 ) )
166	mask_east = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_static_0(k,j,i+1), 0 ) )
167	mask_south = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_static_0(k,j-1,i), 0 ) )
168	mask_north = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_static_0(k,j+1,i), 0 ) )
169
170	tend(k,j,i) = tend(k,j,i) &
171	+ 0.5_wp * ( &
172	mask_east * ( kh(k,j,i) + kh(k,j,i+1) ) &
173	* ( s(k,j,i+1) - s(k,j,i) ) &
174	- mask_west * ( kh(k,j,i) + kh(k,j,i-1) ) &
175	* ( s(k,j,i) - s(k,j,i-1) ) &
176	) * ddx2 * flag &
177	+ 0.5_wp * ( &
178	mask_north * ( kh(k,j,i) + kh(k,j+1,i) ) &
179	* ( s(k,j+1,i) - s(k,j,i) ) &
180	- mask_south * ( kh(k,j,i) + kh(k,j-1,i) ) &
181	* ( s(k,j,i) - s(k,j-1,i) ) &
182	) * ddy2 * flag
183	ENDDO
184
185	!
186	!-- Apply prescribed horizontal wall heatflux where necessary. First,
187	!-- determine start and end index for respective (j,i)-index. Please
188	!-- note, in the flat case following loop will not be entered, as
189	!-- surf_s=1 and surf_e=0. Furtermore, note, no vertical natural surfaces
190	!-- so far.
191	!-- First, for default-type surfaces
192	!-- North-facing vertical default-type surfaces
193	surf_s = surf_def_v(0)%start_index(j,i)
194	surf_e = surf_def_v(0)%end_index(j,i)
195	DO m = surf_s, surf_e
196	k = surf_def_v(0)%k(m)
197	tend(k,j,i) = tend(k,j,i) + s_flux_def_v_north(m) * ddy
198	ENDDO
199	!
200	!-- South-facing vertical default-type surfaces
201	surf_s = surf_def_v(1)%start_index(j,i)
202	surf_e = surf_def_v(1)%end_index(j,i)
203	DO m = surf_s, surf_e
204	k = surf_def_v(1)%k(m)
205	tend(k,j,i) = tend(k,j,i) + s_flux_def_v_south(m) * ddy
206	ENDDO
207	!
208	!-- East-facing vertical default-type surfaces
209	surf_s = surf_def_v(2)%start_index(j,i)
210	surf_e = surf_def_v(2)%end_index(j,i)
211	DO m = surf_s, surf_e
212	k = surf_def_v(2)%k(m)
213	tend(k,j,i) = tend(k,j,i) + s_flux_def_v_east(m) * ddx
214	ENDDO
215	!
216	!-- West-facing vertical default-type surfaces
217	surf_s = surf_def_v(3)%start_index(j,i)
218	surf_e = surf_def_v(3)%end_index(j,i)
219	DO m = surf_s, surf_e
220	k = surf_def_v(3)%k(m)
221	tend(k,j,i) = tend(k,j,i) + s_flux_def_v_west(m) * ddx
222	ENDDO
223	!
224	!-- Now, for natural-type surfaces.
225	!-- North-facing
226	surf_s = surf_lsm_v(0)%start_index(j,i)
227	surf_e = surf_lsm_v(0)%end_index(j,i)
228	DO m = surf_s, surf_e
229	k = surf_lsm_v(0)%k(m)
230	tend(k,j,i) = tend(k,j,i) + s_flux_lsm_v_north(m) * ddy
231	ENDDO
232	!
233	!-- South-facing
234	surf_s = surf_lsm_v(1)%start_index(j,i)
235	surf_e = surf_lsm_v(1)%end_index(j,i)
236	DO m = surf_s, surf_e
237	k = surf_lsm_v(1)%k(m)
238	tend(k,j,i) = tend(k,j,i) + s_flux_lsm_v_south(m) * ddy
239	ENDDO
240	!
241	!-- East-facing
242	surf_s = surf_lsm_v(2)%start_index(j,i)
243	surf_e = surf_lsm_v(2)%end_index(j,i)
244	DO m = surf_s, surf_e
245	k = surf_lsm_v(2)%k(m)
246	tend(k,j,i) = tend(k,j,i) + s_flux_lsm_v_east(m) * ddx
247	ENDDO
248	!
249	!-- West-facing
250	surf_s = surf_lsm_v(3)%start_index(j,i)
251	surf_e = surf_lsm_v(3)%end_index(j,i)
252	DO m = surf_s, surf_e
253	k = surf_lsm_v(3)%k(m)
254	tend(k,j,i) = tend(k,j,i) + s_flux_lsm_v_west(m) * ddx
255	ENDDO
256	!
257	!-- Now, for urban-type surfaces.
258	!-- North-facing
259	surf_s = surf_usm_v(0)%start_index(j,i)
260	surf_e = surf_usm_v(0)%end_index(j,i)
261	DO m = surf_s, surf_e
262	k = surf_usm_v(0)%k(m)
263	tend(k,j,i) = tend(k,j,i) + s_flux_usm_v_north(m) * ddy
264	ENDDO
265	!
266	!-- South-facing
267	surf_s = surf_usm_v(1)%start_index(j,i)
268	surf_e = surf_usm_v(1)%end_index(j,i)
269	DO m = surf_s, surf_e
270	k = surf_usm_v(1)%k(m)
271	tend(k,j,i) = tend(k,j,i) + s_flux_usm_v_south(m) * ddy
272	ENDDO
273	!
274	!-- East-facing
275	surf_s = surf_usm_v(2)%start_index(j,i)
276	surf_e = surf_usm_v(2)%end_index(j,i)
277	DO m = surf_s, surf_e
278	k = surf_usm_v(2)%k(m)
279	tend(k,j,i) = tend(k,j,i) + s_flux_usm_v_east(m) * ddx
280	ENDDO
281	!
282	!-- West-facing
283	surf_s = surf_usm_v(3)%start_index(j,i)
284	surf_e = surf_usm_v(3)%end_index(j,i)
285	DO m = surf_s, surf_e
286	k = surf_usm_v(3)%k(m)
287	tend(k,j,i) = tend(k,j,i) + s_flux_usm_v_west(m) * ddx
288	ENDDO
289
290	!
291	!-- Compute vertical diffusion. In case that surface fluxes have been
292	!-- prescribed or computed at bottom and/or top, index k starts/ends at
293	!-- nzb+2 or nzt-1, respectively. Model top is also mask if top flux
294	!-- is given.
295	DO k = nzb+1, nzt
296	!
297	!-- Determine flags to mask topography below and above. Flag 0 is
298	!-- used to mask topography in general, and flag 8 implies
299	!-- information about use_surface_fluxes. Flag 9 is used to control
300	!-- flux at model top.
301	mask_bottom = MERGE( 1.0_wp, 0.0_wp, &
302	BTEST( wall_flags_static_0(k-1,j,i), 8 ) )
303	mask_top = MERGE( 1.0_wp, 0.0_wp, &
304	BTEST( wall_flags_static_0(k+1,j,i), 8 ) ) * &
305	MERGE( 1.0_wp, 0.0_wp, &
306	BTEST( wall_flags_static_0(k+1,j,i), 9 ) )
307	flag = MERGE( 1.0_wp, 0.0_wp, &
308	BTEST( wall_flags_static_0(k,j,i), 0 ) )
309
310	tend(k,j,i) = tend(k,j,i) &
311	+ 0.5_wp * ( &
312	( kh(k,j,i) + kh(k+1,j,i) ) * &
313	( s(k+1,j,i)-s(k,j,i) ) * ddzu(k+1) &
314	* rho_air_zw(k) &
315	* mask_top &
316	- ( kh(k,j,i) + kh(k-1,j,i) ) * &
317	( s(k,j,i)-s(k-1,j,i) ) * ddzu(k) &
318	* rho_air_zw(k-1) &
319	* mask_bottom &
320	) * ddzw(k) * drho_air(k) &
321	* flag
322	ENDDO
323
324	!
325	!-- Vertical diffusion at horizontal walls.
326	IF ( use_surface_fluxes ) THEN
327	!
328	!-- Default-type surfaces, upward-facing
329	surf_s = surf_def_h(0)%start_index(j,i)
330	surf_e = surf_def_h(0)%end_index(j,i)
331	DO m = surf_s, surf_e
332
333	k = surf_def_h(0)%k(m)
334	tend(k,j,i) = tend(k,j,i) + s_flux_def_h_up(m) &
335	* ddzw(k) * drho_air(k)
336
337	ENDDO
338	!
339	!-- Default-type surfaces, downward-facing
340	surf_s = surf_def_h(1)%start_index(j,i)
341	surf_e = surf_def_h(1)%end_index(j,i)
342	DO m = surf_s, surf_e
343
344	k = surf_def_h(1)%k(m)
345	tend(k,j,i) = tend(k,j,i) + s_flux_def_h_down(m) &
346	* ddzw(k) * drho_air(k)
347
348	ENDDO
349	!
350	!-- Natural-type surfaces, upward-facing
351	surf_s = surf_lsm_h%start_index(j,i)
352	surf_e = surf_lsm_h%end_index(j,i)
353	DO m = surf_s, surf_e
354
355	k = surf_lsm_h%k(m)
356	tend(k,j,i) = tend(k,j,i) + s_flux_lsm_h_up(m) &
357	* ddzw(k) * drho_air(k)
358
359	ENDDO
360	!
361	!-- Urban-type surfaces, upward-facing
362	surf_s = surf_usm_h%start_index(j,i)
363	surf_e = surf_usm_h%end_index(j,i)
364	DO m = surf_s, surf_e
365
366	k = surf_usm_h%k(m)
367	tend(k,j,i) = tend(k,j,i) + s_flux_usm_h_up(m) &
368	* ddzw(k) * drho_air(k)
369
370	ENDDO
371
372	ENDIF
373	!
374	!-- Vertical diffusion at the last computational gridpoint along z-direction
375	IF ( use_top_fluxes ) THEN
376	surf_s = surf_def_h(2)%start_index(j,i)
377	surf_e = surf_def_h(2)%end_index(j,i)
378	DO m = surf_s, surf_e
379
380	k = surf_def_h(2)%k(m)
381	tend(k,j,i) = tend(k,j,i) &
382	+ ( - s_flux_t(m) ) * ddzw(k) * drho_air(k)
383	ENDDO
384	ENDIF
385
386	ENDDO
387	ENDDO
388
389	END SUBROUTINE diffusion_s
390
391	!------------------------------------------------------------------------------!
392	! Description:
393	! ------------
394	!> Call for grid point i,j
395	!------------------------------------------------------------------------------!
396	SUBROUTINE diffusion_s_ij( i, j, s, &
397	s_flux_def_h_up, s_flux_def_h_down, &
398	s_flux_t, &
399	s_flux_lsm_h_up, s_flux_usm_h_up, &
400	s_flux_def_v_north, s_flux_def_v_south, &
401	s_flux_def_v_east, s_flux_def_v_west, &
402	s_flux_lsm_v_north, s_flux_lsm_v_south, &
403	s_flux_lsm_v_east, s_flux_lsm_v_west, &
404	s_flux_usm_v_north, s_flux_usm_v_south, &
405	s_flux_usm_v_east, s_flux_usm_v_west )
406
407	USE arrays_3d, &
408	ONLY: ddzu, ddzw, kh, tend, drho_air, rho_air_zw
409
410	USE control_parameters, &
411	ONLY: use_surface_fluxes, use_top_fluxes
412
413	USE grid_variables, &
414	ONLY: ddx, ddx2, ddy, ddy2
415
416	USE indices, &
417	ONLY: nxlg, nxrg, nyng, nysg, nzb, nzt, wall_flags_static_0
418
419	USE kinds
420
421	USE surface_mod, &
422	ONLY : surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, &
423	surf_usm_v
424
425	IMPLICIT NONE
426
427	INTEGER(iwp) :: i !< running index x direction
428	INTEGER(iwp) :: j !< running index y direction
429	INTEGER(iwp) :: k !< running index z direction
430	INTEGER(iwp) :: m !< running index surface elements
431	INTEGER(iwp) :: surf_e !< End index of surface elements at (j,i)-gridpoint
432	INTEGER(iwp) :: surf_s !< Start index of surface elements at (j,i)-gridpoint
433
434	REAL(wp) :: flag !< flag to mask topography grid points
435	REAL(wp) :: mask_bottom !< flag to mask vertical upward-facing surface
436	REAL(wp) :: mask_east !< flag to mask vertical surface east of the grid point
437	REAL(wp) :: mask_north !< flag to mask vertical surface north of the grid point
438	REAL(wp) :: mask_south !< flag to mask vertical surface south of the grid point
439	REAL(wp) :: mask_west !< flag to mask vertical surface west of the grid point
440	REAL(wp) :: mask_top !< flag to mask vertical downward-facing surface
441
442	REAL(wp), DIMENSION(1:surf_def_v(0)%ns) :: s_flux_def_v_north !< flux at north-facing vertical default-type surfaces
443	REAL(wp), DIMENSION(1:surf_def_v(1)%ns) :: s_flux_def_v_south !< flux at south-facing vertical default-type surfaces
444	REAL(wp), DIMENSION(1:surf_def_v(2)%ns) :: s_flux_def_v_east !< flux at east-facing vertical default-type surfaces
445	REAL(wp), DIMENSION(1:surf_def_v(3)%ns) :: s_flux_def_v_west !< flux at west-facing vertical default-type surfaces
446	REAL(wp), DIMENSION(1:surf_def_h(0)%ns) :: s_flux_def_h_up !< flux at horizontal upward-facing default-type surfaces
447	REAL(wp), DIMENSION(1:surf_def_h(1)%ns) :: s_flux_def_h_down !< flux at horizontal donwward-facing default-type surfaces
448	REAL(wp), DIMENSION(1:surf_lsm_h%ns) :: s_flux_lsm_h_up !< flux at horizontal upward-facing natural-type surfaces
449	REAL(wp), DIMENSION(1:surf_lsm_v(0)%ns) :: s_flux_lsm_v_north !< flux at north-facing vertical urban-type surfaces
450	REAL(wp), DIMENSION(1:surf_lsm_v(1)%ns) :: s_flux_lsm_v_south !< flux at south-facing vertical urban-type surfaces
451	REAL(wp), DIMENSION(1:surf_lsm_v(2)%ns) :: s_flux_lsm_v_east !< flux at east-facing vertical urban-type surfaces
452	REAL(wp), DIMENSION(1:surf_lsm_v(3)%ns) :: s_flux_lsm_v_west !< flux at west-facing vertical urban-type surfaces
453	REAL(wp), DIMENSION(1:surf_usm_h%ns) :: s_flux_usm_h_up !< flux at horizontal upward-facing urban-type surfaces
454	REAL(wp), DIMENSION(1:surf_usm_v(0)%ns) :: s_flux_usm_v_north !< flux at north-facing vertical urban-type surfaces
455	REAL(wp), DIMENSION(1:surf_usm_v(1)%ns) :: s_flux_usm_v_south !< flux at south-facing vertical urban-type surfaces
456	REAL(wp), DIMENSION(1:surf_usm_v(2)%ns) :: s_flux_usm_v_east !< flux at east-facing vertical urban-type surfaces
457	REAL(wp), DIMENSION(1:surf_usm_v(3)%ns) :: s_flux_usm_v_west !< flux at west-facing vertical urban-type surfaces
458	REAL(wp), DIMENSION(1:surf_def_h(2)%ns) :: s_flux_t !< flux at model top
459
460	REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: s !< treated scalar
461
462	!
463	!-- Compute horizontal diffusion
464	DO k = nzb+1, nzt
465	!
466	!-- Predetermine flag to mask topography and wall-bounded grid points
467	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_static_0(k,j,i), 0 ) )
468	!
469	!-- Predetermine flag to mask wall-bounded grid points, equivalent to
470	!-- former s_outer array
471	mask_west = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_static_0(k,j,i-1), 0 ) )
472	mask_east = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_static_0(k,j,i+1), 0 ) )
473	mask_south = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_static_0(k,j-1,i), 0 ) )
474	mask_north = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_static_0(k,j+1,i), 0 ) )
475	!
476	!-- Finally, determine flag to mask both topography itself as well
477	!-- as wall-bounded grid points, which will be treated further below
478
479	tend(k,j,i) = tend(k,j,i) &
480	+ 0.5_wp * ( &
481	mask_east * ( kh(k,j,i) + kh(k,j,i+1) ) &
482	* ( s(k,j,i+1) - s(k,j,i) ) &
483	- mask_west * ( kh(k,j,i) + kh(k,j,i-1) ) &
484	* ( s(k,j,i) - s(k,j,i-1) ) &
485	) * ddx2 * flag &
486	+ 0.5_wp * ( &
487	mask_north * ( kh(k,j,i) + kh(k,j+1,i) ) &
488	* ( s(k,j+1,i) - s(k,j,i) ) &
489	- mask_south * ( kh(k,j,i) + kh(k,j-1,i) ) &
490	* ( s(k,j,i) - s(k,j-1,i) ) &
491	) * ddy2 * flag
492	ENDDO
493
494	!
495	!-- Apply prescribed horizontal wall heatflux where necessary. First,
496	!-- determine start and end index for respective (j,i)-index. Please
497	!-- note, in the flat case following loops will not be entered, as
498	!-- surf_s=1 and surf_e=0. Furtermore, note, no vertical natural surfaces
499	!-- so far.
500	!-- First, for default-type surfaces
501	!-- North-facing vertical default-type surfaces
502	surf_s = surf_def_v(0)%start_index(j,i)
503	surf_e = surf_def_v(0)%end_index(j,i)
504	DO m = surf_s, surf_e
505	k = surf_def_v(0)%k(m)
506	tend(k,j,i) = tend(k,j,i) + s_flux_def_v_north(m) * ddy
507	ENDDO
508	!
509	!-- South-facing vertical default-type surfaces
510	surf_s = surf_def_v(1)%start_index(j,i)
511	surf_e = surf_def_v(1)%end_index(j,i)
512	DO m = surf_s, surf_e
513	k = surf_def_v(1)%k(m)
514	tend(k,j,i) = tend(k,j,i) + s_flux_def_v_south(m) * ddy
515	ENDDO
516	!
517	!-- East-facing vertical default-type surfaces
518	surf_s = surf_def_v(2)%start_index(j,i)
519	surf_e = surf_def_v(2)%end_index(j,i)
520	DO m = surf_s, surf_e
521	k = surf_def_v(2)%k(m)
522	tend(k,j,i) = tend(k,j,i) + s_flux_def_v_east(m) * ddx
523	ENDDO
524	!
525	!-- West-facing vertical default-type surfaces
526	surf_s = surf_def_v(3)%start_index(j,i)
527	surf_e = surf_def_v(3)%end_index(j,i)
528	DO m = surf_s, surf_e
529	k = surf_def_v(3)%k(m)
530	tend(k,j,i) = tend(k,j,i) + s_flux_def_v_west(m) * ddx
531	ENDDO
532	!
533	!-- Now, for natural-type surfaces
534	!-- North-facing
535	surf_s = surf_lsm_v(0)%start_index(j,i)
536	surf_e = surf_lsm_v(0)%end_index(j,i)
537	DO m = surf_s, surf_e
538	k = surf_lsm_v(0)%k(m)
539	tend(k,j,i) = tend(k,j,i) + s_flux_lsm_v_north(m) * ddy
540	ENDDO
541	!
542	!-- South-facing
543	surf_s = surf_lsm_v(1)%start_index(j,i)
544	surf_e = surf_lsm_v(1)%end_index(j,i)
545	DO m = surf_s, surf_e
546	k = surf_lsm_v(1)%k(m)
547	tend(k,j,i) = tend(k,j,i) + s_flux_lsm_v_south(m) * ddy
548	ENDDO
549	!
550	!-- East-facing
551	surf_s = surf_lsm_v(2)%start_index(j,i)
552	surf_e = surf_lsm_v(2)%end_index(j,i)
553	DO m = surf_s, surf_e
554	k = surf_lsm_v(2)%k(m)
555	tend(k,j,i) = tend(k,j,i) + s_flux_lsm_v_east(m) * ddx
556	ENDDO
557	!
558	!-- West-facing
559	surf_s = surf_lsm_v(3)%start_index(j,i)
560	surf_e = surf_lsm_v(3)%end_index(j,i)
561	DO m = surf_s, surf_e
562	k = surf_lsm_v(3)%k(m)
563	tend(k,j,i) = tend(k,j,i) + s_flux_lsm_v_west(m) * ddx
564	ENDDO
565	!
566	!-- Now, for urban-type surfaces
567	!-- North-facing
568	surf_s = surf_usm_v(0)%start_index(j,i)
569	surf_e = surf_usm_v(0)%end_index(j,i)
570	DO m = surf_s, surf_e
571	k = surf_usm_v(0)%k(m)
572	tend(k,j,i) = tend(k,j,i) + s_flux_usm_v_north(m) * ddy
573	ENDDO
574	!
575	!-- South-facing
576	surf_s = surf_usm_v(1)%start_index(j,i)
577	surf_e = surf_usm_v(1)%end_index(j,i)
578	DO m = surf_s, surf_e
579	k = surf_usm_v(1)%k(m)
580	tend(k,j,i) = tend(k,j,i) + s_flux_usm_v_south(m) * ddy
581	ENDDO
582	!
583	!-- East-facing
584	surf_s = surf_usm_v(2)%start_index(j,i)
585	surf_e = surf_usm_v(2)%end_index(j,i)
586	DO m = surf_s, surf_e
587	k = surf_usm_v(2)%k(m)
588	tend(k,j,i) = tend(k,j,i) + s_flux_usm_v_east(m) * ddx
589	ENDDO
590	!
591	!-- West-facing
592	surf_s = surf_usm_v(3)%start_index(j,i)
593	surf_e = surf_usm_v(3)%end_index(j,i)
594	DO m = surf_s, surf_e
595	k = surf_usm_v(3)%k(m)
596	tend(k,j,i) = tend(k,j,i) + s_flux_usm_v_west(m) * ddx
597	ENDDO
598
599
600	!
601	!-- Compute vertical diffusion. In case that surface fluxes have been
602	!-- prescribed or computed at bottom and/or top, index k starts/ends at
603	!-- nzb+2 or nzt-1, respectively. Model top is also mask if top flux
604	!-- is given.
605	DO k = nzb+1, nzt
606	!
607	!-- Determine flags to mask topography below and above. Flag 0 is
608	!-- used to mask topography in general, and flag 8 implies
609	!-- information about use_surface_fluxes. Flag 9 is used to control
610	!-- flux at model top.
611	mask_bottom = MERGE( 1.0_wp, 0.0_wp, &
612	BTEST( wall_flags_static_0(k-1,j,i), 8 ) )
613	mask_top = MERGE( 1.0_wp, 0.0_wp, &
614	BTEST( wall_flags_static_0(k+1,j,i), 8 ) ) * &
615	MERGE( 1.0_wp, 0.0_wp, &
616	BTEST( wall_flags_static_0(k+1,j,i), 9 ) )
617	flag = MERGE( 1.0_wp, 0.0_wp, &
618	BTEST( wall_flags_static_0(k,j,i), 0 ) )
619
620	tend(k,j,i) = tend(k,j,i) &
621	+ 0.5_wp * ( &
622	( kh(k,j,i) + kh(k+1,j,i) ) * &
623	( s(k+1,j,i)-s(k,j,i) ) * ddzu(k+1) &
624	* rho_air_zw(k) &
625	* mask_top &
626	- ( kh(k,j,i) + kh(k-1,j,i) ) * &
627	( s(k,j,i)-s(k-1,j,i) ) * ddzu(k) &
628	* rho_air_zw(k-1) &
629	* mask_bottom &
630	) * ddzw(k) * drho_air(k) &
631	* flag
632	ENDDO
633
634	!
635	!-- Vertical diffusion at horizontal walls.
636	!-- TO DO: Adjust for downward facing walls and mask already in main loop
637	IF ( use_surface_fluxes ) THEN
638	!
639	!-- Default-type surfaces, upward-facing
640	surf_s = surf_def_h(0)%start_index(j,i)
641	surf_e = surf_def_h(0)%end_index(j,i)
642	DO m = surf_s, surf_e
643
644	k = surf_def_h(0)%k(m)
645
646	tend(k,j,i) = tend(k,j,i) + s_flux_def_h_up(m) &
647	* ddzw(k) * drho_air(k)
648	ENDDO
649	!
650	!-- Default-type surfaces, downward-facing
651	surf_s = surf_def_h(1)%start_index(j,i)
652	surf_e = surf_def_h(1)%end_index(j,i)
653	DO m = surf_s, surf_e
654
655	k = surf_def_h(1)%k(m)
656
657	tend(k,j,i) = tend(k,j,i) + s_flux_def_h_down(m) &
658	* ddzw(k) * drho_air(k)
659	ENDDO
660	!
661	!-- Natural-type surfaces, upward-facing
662	surf_s = surf_lsm_h%start_index(j,i)
663	surf_e = surf_lsm_h%end_index(j,i)
664	DO m = surf_s, surf_e
665	k = surf_lsm_h%k(m)
666
667	tend(k,j,i) = tend(k,j,i) + s_flux_lsm_h_up(m) &
668	* ddzw(k) * drho_air(k)
669	ENDDO
670	!
671	!-- Urban-type surfaces, upward-facing
672	surf_s = surf_usm_h%start_index(j,i)
673	surf_e = surf_usm_h%end_index(j,i)
674	DO m = surf_s, surf_e
675	k = surf_usm_h%k(m)
676
677	tend(k,j,i) = tend(k,j,i) + s_flux_usm_h_up(m) &
678	* ddzw(k) * drho_air(k)
679	ENDDO
680	ENDIF
681	!
682	!-- Vertical diffusion at the last computational gridpoint along z-direction
683	IF ( use_top_fluxes ) THEN
684	surf_s = surf_def_h(2)%start_index(j,i)
685	surf_e = surf_def_h(2)%end_index(j,i)
686	DO m = surf_s, surf_e
687
688	k = surf_def_h(2)%k(m)
689	tend(k,j,i) = tend(k,j,i) &
690	+ ( - s_flux_t(m) ) * ddzw(k) * drho_air(k)
691	ENDDO
692	ENDIF
693
694	END SUBROUTINE diffusion_s_ij
695
696	END MODULE diffusion_s_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |