Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |