Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |