Home

Context Navigation

source: palm/trunk/SOURCE/surface_layer_fluxes_mod.f90 @ 3947

Last change on this file since 3947 was 3885, checked in by kanani, 6 years ago
restructure/add location/debug messages
Property svn:keywords set to `Id`
File size: 75.8 KB

Line
1	!> @file surface_layer_fluxes_mod.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	!
21	! Current revisions:
22	! ------------------
23	!
24	!
25	! Former revisions:
26	! -----------------
27	! $Id: surface_layer_fluxes_mod.f90 3885 2019-04-11 11:29:34Z hellstea $
28	! Changes related to global restructuring of location messages and introduction
29	! of additional debug messages
30	!
31	! 3881 2019-04-10 09:31:22Z suehring
32	! Assure that Obukhov length does not become zero
33	!
34	! 3834 2019-03-28 15:40:15Z forkel
35	! added USE chem_gasphase_mod
36	!
37	! 3787 2019-03-07 08:43:54Z raasch
38	! unused variables removed
39	!
40	! 3745 2019-02-15 18:57:56Z suehring
41	! Bugfix, missing calculation of 10cm temperature at vertical building walls,
42	! required for indoor model
43	!
44	! 3744 2019-02-15 18:38:58Z suehring
45	! Some interface calls moved to module_interface + cleanup
46	!
47	! 3668 2019-01-14 12:49:24Z maronga
48	! Removed methods "circular" and "lookup"
49	!
50	! 3655 2019-01-07 16:51:22Z knoop
51	! OpenACC port for SPEC
52	!
53	! 3597 2018-12-04 08:40:18Z maronga
54	! Added routine for calculating near surface air potential temperature (moved
55	! from urban_surface_mod)
56	!
57	! 3547 2018-11-21 13:21:24Z suehring
58	! variables documented
59	!
60	! 3361 2018-10-16 20:39:37Z knoop
61	! Modularization of all bulk cloud physics code components
62	!
63	! 3271 2018-09-24 08:20:34Z suehring
64	! Comment revised
65	!
66	! 3157 2018-07-19 21:08:49Z maronga
67	! Added local free convection velocity scale w_lfc in calculation of uvw_abs.
68	! This can be switche on/off by the user via the flag namelist parameter
69	! use_free_convection_scaling
70	!
71	! 3152 2018-07-19 13:26:52Z suehring
72	! q_surface is now part of surface structure
73	!
74	! 3149 2018-07-19 05:48:09Z maronga
75	! Correct accidental last commit
76	!
77	! 3148 2018-07-19 05:45:25Z maronga
78	! Major bugfix in calculation of Obukhov length
79	!
80	! 3146 2018-07-18 22:36:19Z maronga
81	! move phi_m from turbulence_closure_mod
82	!
83	! 3045 2018-05-28 07:55:41Z Giersch
84	! Error message revised
85	!
86	! 2766 2018-01-22 17:17:47Z kanani
87	! Removed preprocessor directive __chem
88	!
89	! 2718 2018-01-02 08:49:38Z maronga
90	! Corrected "Former revisions" section
91	!
92	! 2696 2017-12-14 17:12:51Z kanani
93	! - Change in file header (GPL part)
94	! - Implementation of chemistry module (FK)
95	! - Added calculation of pt1 and qv1 for all surface types. Added calculation of
96	! pt_surface for default-type surfaces (BM)
97	! - Add flag to disable computation of qsws in case of urban surface (MS)
98	!
99	! 2547 2017-10-16 12:41:56Z schwenkel
100	! extended by cloud_droplets option
101	!
102	! 2321 2017-07-24 15:57:07Z schwenkel
103	! Bugfix: Correct index in lookup table for Obukhov length
104	!
105	! 2299 2017-06-29 10:14:38Z suehring
106	! Adjusted for allow separate spinups of LSM and atmosphere code
107	!
108	! 2292 2017-06-20 09:51:42Z schwenkel
109	! Implementation of new microphysic scheme: cloud_scheme = 'morrison'
110	! includes two more prognostic equations for cloud drop concentration (nc)
111	! and cloud water content (qc).
112	!
113	! 2281 2017-06-13 11:34:50Z suehring
114	! Clean-up unnecessary index access to surface type
115	!
116	! 2233 2017-05-30 18:08:54Z suehring
117	!
118	! 2232 2017-05-30 17:47:52Z suehring
119	! Adjustments to new surface concept
120	! OpenMP bugfix
121	!
122	! 2118 2017-01-17 16:38:49Z raasch
123	! OpenACC directives and related code removed
124	!
125	! 2091 2016-12-21 16:38:18Z suehring
126	! Bugfix in calculation of vsws ( incorrect linear interpolation of us )
127	!
128	! 2076 2016-12-02 13:54:20Z raasch
129	! further openmp bugfix for lookup method
130	!
131	! 2073 2016-11-30 14:34:05Z raasch
132	! openmp bugfix for lookup method
133	!
134	! 2037 2016-10-26 11:15:40Z knoop
135	! Anelastic approximation implemented
136	!
137	! 2011 2016-09-19 17:29:57Z kanani
138	! Flag urban_surface is now defined in module control_parameters.
139	!
140	! 2007 2016-08-24 15:47:17Z kanani
141	! Account for urban surface model in computation of vertical kinematic heatflux
142	!
143	! 2000 2016-08-20 18:09:15Z knoop
144	! Forced header and separation lines into 80 columns
145	!
146	! 1992 2016-08-12 15:14:59Z suehring
147	! Minor bug, declaration of look-up index as INTEGER
148	!
149	! 1960 2016-07-12 16:34:24Z suehring
150	! Treat humidity and passive scalar separately
151	!
152	! 1929 2016-06-09 16:25:25Z suehring
153	! Bugfix: avoid segmentation fault in case one grid point is horizontally
154	! completely surrounded by topography
155	!
156	! 1920 2016-05-30 10:50:15Z suehring
157	! Avoid segmentation fault (see change in 1915) by different initialization of
158	! us instead of adding a very small number in the denominator
159	!
160	! 1915 2016-05-27 11:05:02Z suehring
161	! Bugfix: avoid segmentation fault in case of most_method = 'circular' at first
162	! timestep
163	!
164	! 1850 2016-04-08 13:29:27Z maronga
165	! Module renamed
166	!
167	!
168	! 1822 2016-04-07 07:49:42Z hoffmann
169	! icloud_scheme replaced by microphysics_*
170	!
171	! 1788 2016-03-10 11:01:04Z maronga
172	! Added parameter z0q which replaces z0h in the similarity functions for
173	! humidity.
174	! Syntax layout improved.
175	!
176	! 1757 2016-02-22 15:49:32Z maronga
177	! Minor fixes.
178	!
179	! 1749 2016-02-09 12:19:56Z raasch
180	! further OpenACC adjustments
181	!
182	! 1747 2016-02-08 12:25:53Z raasch
183	! adjustments for OpenACC usage
184	!
185	! 1709 2015-11-04 14:47:01Z maronga
186	! Bugfix: division by zero could occur when calculating rib at low wind speeds
187	! Bugfix: calculation of uv_total for neutral = .T., initial value for ol for
188	! neutral = .T.
189	!
190	! 1705 2015-11-02 14:28:56Z maronga
191	! Typo removed
192	!
193	! 1697 2015-10-28 17:14:10Z raasch
194	! FORTRAN and OpenMP errors removed
195	!
196	! 1696 2015-10-27 10:03:34Z maronga
197	! Modularized and completely re-written version of prandtl_fluxes.f90. In the
198	! course of the re-writing two additional methods have been implemented. See
199	! updated description.
200	!
201	! 1551 2015-03-03 14:18:16Z maronga
202	! Removed land surface model part. The surface fluxes are now always calculated
203	! within prandtl_fluxes, based on the given surface temperature/humidity (which
204	! is either provided by the land surface model, by large scale forcing data, or
205	! directly prescribed by the user.
206	!
207	! 1496 2014-12-02 17:25:50Z maronga
208	! Adapted for land surface model
209	!
210	! 1494 2014-11-21 17:14:03Z maronga
211	! Bugfixes: qs is now calculated before calculation of Rif. calculation of
212	! buoyancy flux in Rif corrected (added missing humidity term), allow use of
213	! topography for coupled runs (not tested)
214	!
215	! 1361 2014-04-16 15:17:48Z hoffmann
216	! Bugfix: calculation of turbulent fluxes of rain water content (qrsws) and rain
217	! drop concentration (nrsws) added
218	!
219	! 1340 2014-03-25 19:45:13Z kanani
220	! REAL constants defined as wp-kind
221	!
222	! 1320 2014-03-20 08:40:49Z raasch
223	! ONLY-attribute added to USE-statements,
224	! kind-parameters added to all INTEGER and REAL declaration statements,
225	! kinds are defined in new module kinds,
226	! old module precision_kind is removed,
227	! revision history before 2012 removed,
228	! comment fields (!:) to be used for variable explanations added to
229	! all variable declaration statements
230	!
231	! 1276 2014-01-15 13:40:41Z heinze
232	! Use LSF_DATA also in case of Dirichlet bottom boundary condition for scalars
233	!
234	! 1257 2013-11-08 15:18:40Z raasch
235	! openACC "kernels do" replaced by "kernels loop", "loop independent" added
236	!
237	! 1036 2012-10-22 13:43:42Z raasch
238	! code put under GPL (PALM 3.9)
239	!
240	! 1015 2012-09-27 09:23:24Z raasch
241	! OpenACC statements added
242	!
243	! 978 2012-08-09 08:28:32Z fricke
244	! roughness length for scalar quantities z0h added
245	!
246	! Revision 1.1 1998/01/23 10:06:06 raasch
247	! Initial revision
248	!
249	!
250	! Description:
251	! ------------
252	!> Diagnostic computation of vertical fluxes in the constant flux layer from the
253	!> values of the variables at grid point k=1 based on Newton iteration
254	!>
255	!> @todo (re)move large_scale_forcing actions
256	!> @todo check/optimize OpenMP directives
257	!> @todo simplify if conditions (which flux need to be computed in which case)
258	!------------------------------------------------------------------------------!
259	MODULE surface_layer_fluxes_mod
260
261	USE arrays_3d, &
262	ONLY: e, kh, nc, nr, pt, q, ql, qc, qr, s, u, v, vpt, w, zu, zw, &
263	drho_air_zw, rho_air_zw, d_exner
264
265	USE basic_constants_and_equations_mod, &
266	ONLY: g, kappa, lv_d_cp, pi, rd_d_rv
267
268	USE chem_gasphase_mod, &
269	ONLY: nvar
270
271	USE chem_modules, &
272	ONLY: constant_csflux
273
274	USE cpulog
275
276	USE control_parameters, &
277	ONLY: air_chemistry, cloud_droplets, &
278	constant_heatflux, constant_scalarflux, &
279	constant_waterflux, coupling_mode, &
280	debug_output, debug_string, &
281	do_output_at_2m, humidity, &
282	ibc_e_b, ibc_pt_b, indoor_model, initializing_actions, &
283	intermediate_timestep_count, intermediate_timestep_count_max, &
284	land_surface, large_scale_forcing, lsf_surf, message_string, &
285	neutral, passive_scalar, pt_surface, q_surface, &
286	run_coupled, surface_pressure, simulated_time, terminate_run, &
287	time_since_reference_point, urban_surface, &
288	use_free_convection_scaling, zeta_max, zeta_min
289
290	USE grid_variables, &
291	ONLY: dx, dy
292
293	USE indices, &
294	ONLY: nxl, nxr, nys, nyn, nzb
295
296	USE kinds
297
298	USE bulk_cloud_model_mod, &
299	ONLY: bulk_cloud_model, microphysics_morrison, microphysics_seifert
300
301	USE pegrid
302
303	USE land_surface_model_mod, &
304	ONLY: aero_resist_kray, skip_time_do_lsm
305
306	USE surface_mod, &
307	ONLY : surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_type, &
308	surf_usm_h, surf_usm_v
309
310
311	IMPLICIT NONE
312
313	INTEGER(iwp) :: i !< loop index x direction
314	INTEGER(iwp) :: j !< loop index y direction
315	INTEGER(iwp) :: k !< loop index z direction
316	INTEGER(iwp) :: l !< loop index for surf type
317
318	INTEGER(iwp), PARAMETER :: num_steps = 15000 !< number of steps in the lookup table
319
320	LOGICAL :: coupled_run !< Flag for coupled atmosphere-ocean runs
321	LOGICAL :: downward = .FALSE.!< Flag indicating downward-facing horizontal surface
322	LOGICAL :: mom_uv = .FALSE. !< Flag indicating calculation of usvs and vsus at vertical surfaces
323	LOGICAL :: mom_w = .FALSE. !< Flag indicating calculation of wsus and wsvs at vertical surfaces
324	LOGICAL :: mom_tke = .FALSE. !< Flag indicating calculation of momentum fluxes at vertical surfaces used for TKE production
325	LOGICAL :: surf_vertical !< Flag indicating vertical surfaces
326
327	REAL(wp) :: e_s, & !< Saturation water vapor pressure
328	ol_max = 1.0E6_wp, & !< Maximum Obukhov length
329	z_mo !< Height of the constant flux layer where MOST is assumed
330
331	TYPE(surf_type), POINTER :: surf !< surf-type array, used to generalize subroutines
332
333
334	SAVE
335
336	PRIVATE
337
338	PUBLIC init_surface_layer_fluxes, phi_m, surface_layer_fluxes
339
340	INTERFACE init_surface_layer_fluxes
341	MODULE PROCEDURE init_surface_layer_fluxes
342	END INTERFACE init_surface_layer_fluxes
343
344	INTERFACE phi_m
345	MODULE PROCEDURE phi_m
346	END INTERFACE phi_m
347
348	INTERFACE surface_layer_fluxes
349	MODULE PROCEDURE surface_layer_fluxes
350	END INTERFACE surface_layer_fluxes
351
352
353	CONTAINS
354
355
356	!------------------------------------------------------------------------------!
357	! Description:
358	! ------------
359	!> Main routine to compute the surface fluxes
360	!------------------------------------------------------------------------------!
361	SUBROUTINE surface_layer_fluxes
362
363	IMPLICIT NONE
364
365
366	IF ( debug_output ) CALL debug_message( 'surface_layer_fluxes', 'start' )
367
368	surf_vertical = .FALSE. !< flag indicating vertically orientated surface elements
369	downward = .FALSE. !< flag indicating downward-facing surface elements
370	!
371	!-- Derive potential temperature and specific humidity at first grid level
372	!-- from the fields pt and q
373	!
374	!-- First call for horizontal default-type surfaces (l=0 - upward facing,
375	!-- l=1 - downward facing)
376	DO l = 0, 1
377	IF ( surf_def_h(l)%ns >= 1 ) THEN
378	surf => surf_def_h(l)
379	CALL calc_pt_q
380	IF ( .NOT. neutral ) THEN
381	CALL calc_pt_surface
382	IF ( humidity ) THEN
383	CALL calc_q_surface
384	CALL calc_vpt_surface
385	ENDIF
386	ENDIF
387	ENDIF
388	ENDDO
389	!
390	!-- Call for natural-type horizontal surfaces
391	IF ( surf_lsm_h%ns >= 1 ) THEN
392	surf => surf_lsm_h
393	CALL calc_pt_q
394	ENDIF
395
396	!
397	!-- Call for urban-type horizontal surfaces
398	IF ( surf_usm_h%ns >= 1 ) THEN
399	surf => surf_usm_h
400	CALL calc_pt_q
401	ENDIF
402
403	!
404	!-- Call for natural-type vertical surfaces
405	DO l = 0, 3
406	IF ( surf_lsm_v(l)%ns >= 1 ) THEN
407	surf => surf_lsm_v(l)
408	CALL calc_pt_q
409	ENDIF
410
411	!-- Call for urban-type vertical surfaces
412	IF ( surf_usm_v(l)%ns >= 1 ) THEN
413	surf => surf_usm_v(l)
414	CALL calc_pt_q
415	ENDIF
416	ENDDO
417
418	!
419	!-- First, calculate the new Obukhov length, then new friction velocity,
420	!-- followed by the new scaling parameters (th, q, etc.), and the new
421	!-- surface fluxes if required. Note, each routine is called for different surface types.
422	!-- First call for default-type horizontal surfaces, for natural- and
423	!-- urban-type surfaces. Note, at this place only upward-facing horizontal
424	!-- surfaces are treated.
425
426	!
427	!-- Default-type upward-facing horizontal surfaces
428	IF ( surf_def_h(0)%ns >= 1 ) THEN
429	surf => surf_def_h(0)
430	CALL calc_uvw_abs
431	IF ( .NOT. neutral ) CALL calc_ol
432	CALL calc_us
433	CALL calc_scaling_parameters
434	CALL calc_surface_fluxes
435	IF ( do_output_at_2m ) THEN
436	CALL calc_pt_near_surface ( '2m' )
437	ENDIF
438	ENDIF
439	!
440	!-- Natural-type horizontal surfaces
441	IF ( surf_lsm_h%ns >= 1 ) THEN
442	surf => surf_lsm_h
443	CALL calc_uvw_abs
444	IF ( .NOT. neutral ) CALL calc_ol
445	CALL calc_us
446	CALL calc_scaling_parameters
447	CALL calc_surface_fluxes
448	IF ( do_output_at_2m ) THEN
449	CALL calc_pt_near_surface ( '2m' )
450	ENDIF
451	ENDIF
452	!
453	!-- Urban-type horizontal surfaces
454	IF ( surf_usm_h%ns >= 1 ) THEN
455	surf => surf_usm_h
456	CALL calc_uvw_abs
457	IF ( .NOT. neutral ) CALL calc_ol
458	CALL calc_us
459	CALL calc_scaling_parameters
460	CALL calc_surface_fluxes
461	IF ( do_output_at_2m ) THEN
462	CALL calc_pt_near_surface ( '2m' )
463	ENDIF
464	!
465	!-- Calculate 10cm temperature, required in indoor model
466	IF ( indoor_model ) CALL calc_pt_near_surface ( '10cm' )
467	ENDIF
468
469	!
470	!-- Treat downward-facing horizontal surfaces. Note, so far, these are
471	!-- always default type. Stratification is not considered
472	!-- in this case, hence, no further distinction between different
473	!-- most_method is required.
474	IF ( surf_def_h(1)%ns >= 1 ) THEN
475	downward = .TRUE.
476	surf => surf_def_h(1)
477	CALL calc_uvw_abs
478	CALL calc_us
479	CALL calc_surface_fluxes
480	downward = .FALSE.
481	ENDIF
482	!
483	!-- Calculate surfaces fluxes at vertical surfaces for momentum
484	!-- and subgrid-scale TKE.
485	!-- No stability is considered. Therefore, scaling parameters and Obukhov-
486	!-- length do not need to be calculated and no distinction in 'circular',
487	!-- 'Newton' or 'lookup' is necessary so far.
488	!-- Note, this will change if stability is once considered.
489	surf_vertical = .TRUE.
490	!
491	!-- Calculate horizontal momentum fluxes at north- and south-facing
492	!-- surfaces(usvs).
493	!-- For default-type surfaces
494	mom_uv = .TRUE.
495	DO l = 0, 1
496	IF ( surf_def_v(l)%ns >= 1 ) THEN
497	surf => surf_def_v(l)
498	!
499	!-- Compute surface-parallel velocity
500	CALL calc_uvw_abs_v_ugrid
501	!
502	!-- Compute respective friction velocity on staggered grid
503	CALL calc_us
504	!
505	!-- Compute respective surface fluxes for momentum and TKE
506	CALL calc_surface_fluxes
507	ENDIF
508	ENDDO
509	!
510	!-- For natural-type surfaces. Please note, even though stability is not
511	!-- considered for the calculation of momentum fluxes at vertical surfaces,
512	!-- scaling parameters and Obukov length are calculated nevertheless in this
513	!-- case. This is due to the requirement of ts in parameterization of heat
514	!-- flux in land-surface model in case of aero_resist_kray is not true.
515	IF ( .NOT. aero_resist_kray ) THEN
516	DO l = 0, 1
517	IF ( surf_lsm_v(l)%ns >= 1 ) THEN
518	surf => surf_lsm_v(l)
519	!
520	!-- Compute surface-parallel velocity
521	CALL calc_uvw_abs_v_ugrid
522	!
523	!-- Compute Obukhov length
524	IF ( .NOT. neutral ) CALL calc_ol
525	!
526	!-- Compute respective friction velocity on staggered grid
527	CALL calc_us
528	!
529	!-- Compute scaling parameters
530	CALL calc_scaling_parameters
531	!
532	!-- Compute respective surface fluxes for momentum and TKE
533	CALL calc_surface_fluxes
534	ENDIF
535	ENDDO
536	!
537	!-- No ts is required, so scaling parameters and Obukhov length do not need
538	!-- to be computed.
539	ELSE
540	DO l = 0, 1
541	IF ( surf_lsm_v(l)%ns >= 1 ) THEN
542	surf => surf_lsm_v(l)
543	!
544	!-- Compute surface-parallel velocity
545	CALL calc_uvw_abs_v_ugrid
546	!
547	!-- Compute respective friction velocity on staggered grid
548	CALL calc_us
549	!
550	!-- Compute respective surface fluxes for momentum and TKE
551	CALL calc_surface_fluxes
552	ENDIF
553	ENDDO
554	ENDIF
555	!
556	!-- For urban-type surfaces
557	DO l = 0, 1
558	IF ( surf_usm_v(l)%ns >= 1 ) THEN
559	surf => surf_usm_v(l)
560	!
561	!-- Compute surface-parallel velocity
562	CALL calc_uvw_abs_v_ugrid
563	!
564	!-- Compute respective friction velocity on staggered grid
565	CALL calc_us
566	!
567	!-- Compute respective surface fluxes for momentum and TKE
568	CALL calc_surface_fluxes
569	!
570	!-- Calculate 10cm temperature, required in indoor model
571	IF ( indoor_model ) CALL calc_pt_near_surface ( '10cm' )
572	ENDIF
573	ENDDO
574	!
575	!-- Calculate horizontal momentum fluxes at east- and west-facing
576	!-- surfaces (vsus).
577	!-- For default-type surfaces
578	DO l = 2, 3
579	IF ( surf_def_v(l)%ns >= 1 ) THEN
580	surf => surf_def_v(l)
581	!
582	!-- Compute surface-parallel velocity
583	CALL calc_uvw_abs_v_vgrid
584	!
585	!-- Compute respective friction velocity on staggered grid
586	CALL calc_us
587	!
588	!-- Compute respective surface fluxes for momentum and TKE
589	CALL calc_surface_fluxes
590
591	ENDIF
592	ENDDO
593	!
594	!-- For natural-type surfaces. Please note, even though stability is not
595	!-- considered for the calculation of momentum fluxes at vertical surfaces,
596	!-- scaling parameters and Obukov length are calculated nevertheless in this
597	!-- case. This is due to the requirement of ts in parameterization of heat
598	!-- flux in land-surface model in case of aero_resist_kray is not true.
599	IF ( .NOT. aero_resist_kray ) THEN
600	DO l = 2, 3
601	IF ( surf_lsm_v(l)%ns >= 1 ) THEN
602	surf => surf_lsm_v(l)
603	!
604	!-- Compute surface-parallel velocity
605	CALL calc_uvw_abs_v_vgrid
606	!
607	!-- Compute Obukhov length
608	IF ( .NOT. neutral ) CALL calc_ol
609	!
610	!-- Compute respective friction velocity on staggered grid
611	CALL calc_us
612	!
613	!-- Compute scaling parameters
614	CALL calc_scaling_parameters
615	!
616	!-- Compute respective surface fluxes for momentum and TKE
617	CALL calc_surface_fluxes
618	ENDIF
619	ENDDO
620	ELSE
621	DO l = 2, 3
622	IF ( surf_lsm_v(l)%ns >= 1 ) THEN
623	surf => surf_lsm_v(l)
624	!
625	!-- Compute surface-parallel velocity
626	CALL calc_uvw_abs_v_vgrid
627	!
628	!-- Compute respective friction velocity on staggered grid
629	CALL calc_us
630	!
631	!-- Compute respective surface fluxes for momentum and TKE
632	CALL calc_surface_fluxes
633	ENDIF
634	ENDDO
635	ENDIF
636	!
637	!-- For urban-type surfaces
638	DO l = 2, 3
639	IF ( surf_usm_v(l)%ns >= 1 ) THEN
640	surf => surf_usm_v(l)
641	!
642	!-- Compute surface-parallel velocity
643	CALL calc_uvw_abs_v_vgrid
644	!
645	!-- Compute respective friction velocity on staggered grid
646	CALL calc_us
647	!
648	!-- Compute respective surface fluxes for momentum and TKE
649	CALL calc_surface_fluxes
650	!
651	!-- Calculate 10cm temperature, required in indoor model
652	IF ( indoor_model ) CALL calc_pt_near_surface ( '10cm' )
653	ENDIF
654	ENDDO
655	mom_uv = .FALSE.
656	!
657	!-- Calculate horizontal momentum fluxes of w (wsus and wsvs) at vertial
658	!-- surfaces.
659	mom_w = .TRUE.
660	!
661	!-- Default-type surfaces
662	DO l = 0, 3
663	IF ( surf_def_v(l)%ns >= 1 ) THEN
664	surf => surf_def_v(l)
665	CALL calc_uvw_abs_v_wgrid
666	CALL calc_us
667	CALL calc_surface_fluxes
668	ENDIF
669	ENDDO
670	!
671	!-- Natural-type surfaces
672	DO l = 0, 3
673	IF ( surf_lsm_v(l)%ns >= 1 ) THEN
674	surf => surf_lsm_v(l)
675	CALL calc_uvw_abs_v_wgrid
676	CALL calc_us
677	CALL calc_surface_fluxes
678	ENDIF
679	ENDDO
680	!
681	!-- Urban-type surfaces
682	DO l = 0, 3
683	IF ( surf_usm_v(l)%ns >= 1 ) THEN
684	surf => surf_usm_v(l)
685	CALL calc_uvw_abs_v_wgrid
686	CALL calc_us
687	CALL calc_surface_fluxes
688	ENDIF
689	ENDDO
690	mom_w = .FALSE.
691	!
692	!-- Calculate momentum fluxes usvs, vsus, wsus and wsvs at vertical
693	!-- surfaces for TKE production. Note, here, momentum fluxes are defined
694	!-- at grid center and are not staggered as before.
695	mom_tke = .TRUE.
696	!
697	!-- Default-type surfaces
698	DO l = 0, 3
699	IF ( surf_def_v(l)%ns >= 1 ) THEN
700	surf => surf_def_v(l)
701	CALL calc_uvw_abs_v_sgrid
702	CALL calc_us
703	CALL calc_surface_fluxes
704	ENDIF
705	ENDDO
706	!
707	!-- Natural-type surfaces
708	DO l = 0, 3
709	IF ( surf_lsm_v(l)%ns >= 1 ) THEN
710	surf => surf_lsm_v(l)
711	CALL calc_uvw_abs_v_sgrid
712	CALL calc_us
713	CALL calc_surface_fluxes
714	ENDIF
715	ENDDO
716	!
717	!-- Urban-type surfaces
718	DO l = 0, 3
719	IF ( surf_usm_v(l)%ns >= 1 ) THEN
720	surf => surf_usm_v(l)
721	CALL calc_uvw_abs_v_sgrid
722	CALL calc_us
723	CALL calc_surface_fluxes
724	ENDIF
725	ENDDO
726	mom_tke = .FALSE.
727
728	IF ( debug_output ) CALL debug_message( 'surface_layer_fluxes', 'end' )
729
730	END SUBROUTINE surface_layer_fluxes
731
732
733	!------------------------------------------------------------------------------!
734	! Description:
735	! ------------
736	!> Initializing actions for the surface layer routine. Basically, this involves
737	!> the preparation of a lookup table for the the bulk Richardson number vs
738	!> Obukhov length L when using the lookup table method.
739	!------------------------------------------------------------------------------!
740	SUBROUTINE init_surface_layer_fluxes
741
742	IMPLICIT NONE
743
744
745	CALL location_message( 'initializing surface layer', 'start' )
746
747	!
748	!-- In case of runs with neutral statification, set Obukhov length to a
749	!-- large value
750	IF ( neutral ) THEN
751	IF ( surf_def_h(0)%ns >= 1 ) surf_def_h(0)%ol = 1.0E10_wp
752	IF ( surf_lsm_h%ns >= 1 ) surf_lsm_h%ol = 1.0E10_wp
753	IF ( surf_usm_h%ns >= 1 ) surf_usm_h%ol = 1.0E10_wp
754	ENDIF
755
756	CALL location_message( 'initializing surface layer', 'finished' )
757
758	END SUBROUTINE init_surface_layer_fluxes
759
760
761	!------------------------------------------------------------------------------!
762	! Description:
763	! ------------
764	!> Compute the absolute value of the horizontal velocity (relative to the
765	!> surface) for horizontal surface elements. This is required by all methods.
766	!------------------------------------------------------------------------------!
767	SUBROUTINE calc_uvw_abs
768
769	IMPLICIT NONE
770
771	INTEGER(iwp) :: i !< running index x direction
772	INTEGER(iwp) :: ibit !< flag to mask computation of relative velocity in case of downward-facing surfaces
773	INTEGER(iwp) :: j !< running index y direction
774	INTEGER(iwp) :: k !< running index z direction
775	INTEGER(iwp) :: m !< running index surface elements
776
777	REAL(wp) :: w_lfc !< local free convection velocity scale
778	!
779	!-- ibit is 1 for upward-facing surfaces, zero for downward-facing surfaces.
780	ibit = MERGE( 1, 0, .NOT. downward )
781
782	!$ACC PARALLEL LOOP PRIVATE(i, j, k, w_lfc) &
783	!$ACC PRESENT(surf, u, v)
784	DO m = 1, surf%ns
785
786	i = surf%i(m)
787	j = surf%j(m)
788	k = surf%k(m)
789
790	!
791	!-- Calculate free convection velocity scale w_lfc is
792	!-- use_free_convection_scaling = .T.. This will maintain a horizontal
793	!-- velocity even for very weak wind convective conditions. SIGN is used
794	!-- to set w_lfc to zero under stable conditions.
795	IF ( use_free_convection_scaling ) THEN
796	w_lfc = ABS(g / surf%pt1(m) * surf%z_mo(m) * surf%shf(m))
797	w_lfc = ( 0.5_wp * ( w_lfc + SIGN(w_lfc,surf%shf(m)) ) )**(0.33333_wp)
798	ELSE
799	w_lfc = 0.0_wp
800	ENDIF
801
802	!
803	!-- Compute the absolute value of the horizontal velocity.
804	!-- (relative to the surface in case the lower surface is the ocean).
805	!-- Please note, in new surface modelling concept the index values changed,
806	!-- i.e. the reference grid point is not the surface-grid point itself but
807	!-- the first grid point outside of the topography.
808	!-- Note, in case of coupled ocean-atmosphere simulations relative velocity
809	!-- with respect to the ocean surface is used, hence, (k-1,j,i) values
810	!-- are used to calculate the absolute velocity.
811	!-- However, this do not apply for downward-facing walls. To mask this,
812	!-- use ibit, which checks for upward/downward-facing surfaces.
813	surf%uvw_abs(m) = SQRT( &
814	( 0.5_wp * ( u(k,j,i) + u(k,j,i+1) &
815	- ( u(k-1,j,i) + u(k-1,j,i+1) &
816	) * ibit &
817	) &
818	)**2 + &
819	( 0.5_wp * ( v(k,j,i) + v(k,j+1,i) &
820	- ( v(k-1,j,i) + v(k-1,j+1,i) &
821	) * ibit &
822	) &
823	)2 + w_lfc2 &
824	)
825
826
827
828	ENDDO
829
830	END SUBROUTINE calc_uvw_abs
831
832
833	!------------------------------------------------------------------------------!
834	! Description:
835	! ------------
836	!> Compute the absolute value of the horizontal velocity (relative to the
837	!> surface) for horizontal surface elements. This is required by all methods.
838	!------------------------------------------------------------------------------!
839	SUBROUTINE calc_uvw_abs_v_ugrid
840
841	IMPLICIT NONE
842
843	INTEGER(iwp) :: i !< running index x direction
844	INTEGER(iwp) :: j !< running index y direction
845	INTEGER(iwp) :: k !< running index z direction
846	INTEGER(iwp) :: m !< running index surface elements
847
848	REAL(wp) :: u_i !< u-component on xu-grid
849	REAL(wp) :: w_i !< w-component on xu-grid
850
851
852	DO m = 1, surf%ns
853	i = surf%i(m)
854	j = surf%j(m)
855	k = surf%k(m)
856	!
857	!-- Compute the absolute value of the surface parallel velocity on u-grid.
858	u_i = u(k,j,i)
859	w_i = 0.25_wp * ( w(k-1,j,i-1) + w(k-1,j,i) + &
860	w(k,j,i-1) + w(k,j,i) )
861
862	surf%uvw_abs(m) = SQRT( u_i2 + w_i2 )
863
864	ENDDO
865
866	END SUBROUTINE calc_uvw_abs_v_ugrid
867
868	!------------------------------------------------------------------------------!
869	! Description:
870	! ------------
871	!> Compute the absolute value of the horizontal velocity (relative to the
872	!> surface) for horizontal surface elements. This is required by all methods.
873	!------------------------------------------------------------------------------!
874	SUBROUTINE calc_uvw_abs_v_vgrid
875
876	IMPLICIT NONE
877
878	INTEGER(iwp) :: i !< running index x direction
879	INTEGER(iwp) :: j !< running index y direction
880	INTEGER(iwp) :: k !< running index z direction
881	INTEGER(iwp) :: m !< running index surface elements
882
883	REAL(wp) :: v_i !< v-component on yv-grid
884	REAL(wp) :: w_i !< w-component on yv-grid
885
886
887	DO m = 1, surf%ns
888	i = surf%i(m)
889	j = surf%j(m)
890	k = surf%k(m)
891
892	v_i = u(k,j,i)
893	w_i = 0.25_wp * ( w(k-1,j-1,i) + w(k-1,j,i) + &
894	w(k,j-1,i) + w(k,j,i) )
895
896	surf%uvw_abs(m) = SQRT( v_i2 + w_i2 )
897
898	ENDDO
899
900	END SUBROUTINE calc_uvw_abs_v_vgrid
901
902	!------------------------------------------------------------------------------!
903	! Description:
904	! ------------
905	!> Compute the absolute value of the horizontal velocity (relative to the
906	!> surface) for horizontal surface elements. This is required by all methods.
907	!------------------------------------------------------------------------------!
908	SUBROUTINE calc_uvw_abs_v_wgrid
909
910	IMPLICIT NONE
911
912	INTEGER(iwp) :: i !< running index x direction
913	INTEGER(iwp) :: j !< running index y direction
914	INTEGER(iwp) :: k !< running index z direction
915	INTEGER(iwp) :: m !< running index surface elements
916
917	REAL(wp) :: u_i !< u-component on x-zw-grid
918	REAL(wp) :: v_i !< v-component on y-zw-grid
919	REAL(wp) :: w_i !< w-component on zw-grid
920	!
921	!-- North- (l=0) and south-facing (l=1) surfaces
922	IF ( l == 0 .OR. l == 1 ) THEN
923	DO m = 1, surf%ns
924	i = surf%i(m)
925	j = surf%j(m)
926	k = surf%k(m)
927
928	u_i = 0.25_wp * ( u(k+1,j,i+1) + u(k+1,j,i) + &
929	u(k,j,i+1) + u(k,j,i) )
930	v_i = 0.0_wp
931	w_i = w(k,j,i)
932
933	surf%uvw_abs(m) = SQRT( u_i2 + v_i2 + w_i**2 )
934	ENDDO
935	!
936	!-- East- (l=2) and west-facing (l=3) surfaces
937	ELSE
938	DO m = 1, surf%ns
939	i = surf%i(m)
940	j = surf%j(m)
941	k = surf%k(m)
942
943	u_i = 0.0_wp
944	v_i = 0.25_wp * ( v(k+1,j+1,i) + v(k+1,j,i) + &
945	v(k,j+1,i) + v(k,j,i) )
946	w_i = w(k,j,i)
947
948	surf%uvw_abs(m) = SQRT( u_i2 + v_i2 + w_i**2 )
949	ENDDO
950	ENDIF
951
952	END SUBROUTINE calc_uvw_abs_v_wgrid
953
954	!------------------------------------------------------------------------------!
955	! Description:
956	! ------------
957	!> Compute the absolute value of the horizontal velocity (relative to the
958	!> surface) for horizontal surface elements. This is required by all methods.
959	!------------------------------------------------------------------------------!
960	SUBROUTINE calc_uvw_abs_v_sgrid
961
962	IMPLICIT NONE
963
964	INTEGER(iwp) :: i !< running index x direction
965	INTEGER(iwp) :: j !< running index y direction
966	INTEGER(iwp) :: k !< running index z direction
967	INTEGER(iwp) :: m !< running index surface elements
968
969	REAL(wp) :: u_i !< u-component on scalar grid
970	REAL(wp) :: v_i !< v-component on scalar grid
971	REAL(wp) :: w_i !< w-component on scalar grid
972
973	!
974	!-- North- (l=0) and south-facing (l=1) walls
975	IF ( l == 0 .OR. l == 1 ) THEN
976	DO m = 1, surf%ns
977	i = surf%i(m)
978	j = surf%j(m)
979	k = surf%k(m)
980
981	u_i = 0.5_wp * ( u(k,j,i) + u(k,j,i+1) )
982	v_i = 0.0_wp
983	w_i = 0.5_wp * ( w(k,j,i) + w(k-1,j,i) )
984
985	surf%uvw_abs(m) = SQRT( u_i2 + v_i2 + w_i**2 )
986	ENDDO
987	!
988	!-- East- (l=2) and west-facing (l=3) walls
989	ELSE
990	DO m = 1, surf%ns
991	i = surf%i(m)
992	j = surf%j(m)
993	k = surf%k(m)
994
995	u_i = 0.0_wp
996	v_i = 0.5_wp * ( v(k,j,i) + v(k,j+1,i) )
997	w_i = 0.5_wp * ( w(k,j,i) + w(k-1,j,i) )
998
999	surf%uvw_abs(m) = SQRT( u_i2 + v_i2 + w_i**2 )
1000	ENDDO
1001	ENDIF
1002
1003	END SUBROUTINE calc_uvw_abs_v_sgrid
1004
1005
1006	!------------------------------------------------------------------------------!
1007	! Description:
1008	! ------------
1009	!> Calculate the Obukhov length (L) and Richardson flux number (z/L)
1010	!------------------------------------------------------------------------------!
1011	SUBROUTINE calc_ol
1012
1013	IMPLICIT NONE
1014
1015	INTEGER(iwp) :: iter !< Newton iteration step
1016	INTEGER(iwp) :: m !< loop variable over all horizontal wall elements
1017
1018	REAL(wp) :: f, & !< Function for Newton iteration: f = Ri - [...]/[...]^2 = 0
1019	f_d_ol, & !< Derivative of f
1020	ol_l, & !< Lower bound of L for Newton iteration
1021	ol_m, & !< Previous value of L for Newton iteration
1022	ol_old, & !< Previous time step value of L
1023	ol_u !< Upper bound of L for Newton iteration
1024
1025	!
1026	!-- Evaluate bulk Richardson number (calculation depends on
1027	!-- definition based on setting of boundary conditions
1028	IF ( ibc_pt_b /= 1 ) THEN
1029	IF ( humidity ) THEN
1030	!$OMP PARALLEL DO PRIVATE( z_mo )
1031	DO m = 1, surf%ns
1032
1033	z_mo = surf%z_mo(m)
1034
1035	surf%rib(m) = g * z_mo &
1036	* ( surf%vpt1(m) - surf%vpt_surface(m) ) &
1037	/ ( surf%uvw_abs(m)*2 surf%vpt1(m) &
1038	+ 1.0E-20_wp )
1039	ENDDO
1040	ELSE
1041	!$OMP PARALLEL DO PRIVATE( z_mo )
1042	DO m = 1, surf%ns
1043
1044	z_mo = surf%z_mo(m)
1045
1046	surf%rib(m) = g * z_mo &
1047	* ( surf%pt1(m) - surf%pt_surface(m) ) &
1048	/ ( surf%uvw_abs(m)*2 surf%pt1(m) + 1.0E-20_wp )
1049	ENDDO
1050	ENDIF
1051	ELSE
1052	IF ( humidity ) THEN
1053	!$OMP PARALLEL DO PRIVATE( k, z_mo )
1054	DO m = 1, surf%ns
1055
1056	k = surf%k(m)
1057
1058	z_mo = surf%z_mo(m)
1059
1060	surf%rib(m) = - g * z_mo * ( ( 1.0_wp + 0.61_wp &
1061	* surf%qv1(m) ) * surf%shf(m) + 0.61_wp &
1062	* surf%pt1(m) * surf%qsws(m) ) * &
1063	drho_air_zw(k-1) / &
1064	( surf%uvw_abs(m)*3 surf%vpt1(m) * kappa**2 &
1065	+ 1.0E-20_wp )
1066	ENDDO
1067	ELSE
1068	!$OMP PARALLEL DO PRIVATE( k, z_mo )
1069	!$ACC PARALLEL LOOP PRIVATE(k, z_mo) &
1070	!$ACC PRESENT(surf, drho_air_zw)
1071	DO m = 1, surf%ns
1072
1073	k = surf%k(m)
1074
1075	z_mo = surf%z_mo(m)
1076
1077	surf%rib(m) = - g * z_mo * surf%shf(m) * &
1078	drho_air_zw(k-1) / &
1079	( surf%uvw_abs(m)*3 surf%pt1(m) * kappa**2 &
1080	+ 1.0E-20_wp )
1081	ENDDO
1082	ENDIF
1083	ENDIF
1084
1085	!
1086	!-- Calculate the Obukhov length using Newton iteration
1087	!$ACC PARALLEL LOOP PRIVATE(i, j, z_mo) &
1088	!$ACC PRIVATE(ol_old, ol_m, ol_l, ol_u, f, f_d_ol) &
1089	!$ACC PRESENT(surf)
1090	DO m = 1, surf%ns
1091
1092	i = surf%i(m)
1093	j = surf%j(m)
1094
1095	z_mo = surf%z_mo(m)
1096
1097	!
1098	!-- Store current value in case the Newton iteration fails
1099	ol_old = surf%ol(m)
1100
1101	!
1102	!-- Ensure that the bulk Richardson number and the Obukhov
1103	!-- length have the same sign
1104	IF ( surf%rib(m) * surf%ol(m) < 0.0_wp .OR. &
1105	ABS( surf%ol(m) ) == ol_max ) THEN
1106	IF ( surf%rib(m) > 1.0_wp ) surf%ol(m) = 0.01_wp
1107	IF ( surf%rib(m) < 0.0_wp ) surf%ol(m) = -0.01_wp
1108	ENDIF
1109	!
1110	!-- Iteration to find Obukhov length
1111	iter = 0
1112	DO
1113	iter = iter + 1
1114	!
1115	!-- In case of divergence, use the value of the previous time step
1116	IF ( iter > 1000 ) THEN
1117	surf%ol(m) = ol_old
1118	EXIT
1119	ENDIF
1120
1121	ol_m = surf%ol(m)
1122	ol_l = ol_m - 0.001_wp * ol_m
1123	ol_u = ol_m + 0.001_wp * ol_m
1124
1125
1126	IF ( ibc_pt_b /= 1 ) THEN
1127	!
1128	!-- Calculate f = Ri - [...]/[...]^2 = 0
1129	f = surf%rib(m) - ( z_mo / ol_m ) * ( &
1130	LOG( z_mo / surf%z0h(m) ) &
1131	- psi_h( z_mo / ol_m ) &
1132	+ psi_h( surf%z0h(m) / &
1133	ol_m ) &
1134	) &
1135	/ ( LOG( z_mo / surf%z0(m) ) &
1136	- psi_m( z_mo / ol_m ) &
1137	+ psi_m( surf%z0(m) / ol_m ) &
1138	)**2
1139
1140	!
1141	!-- Calculate df/dL
1142	f_d_ol = ( - ( z_mo / ol_u ) * ( LOG( z_mo / &
1143	surf%z0h(m) ) &
1144	- psi_h( z_mo / ol_u ) &
1145	+ psi_h( surf%z0h(m) / ol_u ) &
1146	) &
1147	/ ( LOG( z_mo / surf%z0(m) ) &
1148	- psi_m( z_mo / ol_u ) &
1149	+ psi_m( surf%z0(m) / ol_u ) &
1150	)**2 &
1151	+ ( z_mo / ol_l ) * ( LOG( z_mo / surf%z0h(m) ) &
1152	- psi_h( z_mo / ol_l ) &
1153	+ psi_h( surf%z0h(m) / ol_l ) &
1154	) &
1155	/ ( LOG( z_mo / surf%z0(m) ) &
1156	- psi_m( z_mo / ol_l ) &
1157	+ psi_m( surf%z0(m) / ol_l ) &
1158	)**2 &
1159	) / ( ol_u - ol_l )
1160	ELSE
1161	!
1162	!-- Calculate f = Ri - 1 /[...]^3 = 0
1163	f = surf%rib(m) - ( z_mo / ol_m ) / &
1164	( LOG( z_mo / surf%z0(m) ) &
1165	- psi_m( z_mo / ol_m ) &
1166	+ psi_m( surf%z0(m) / ol_m ) &
1167	)**3
1168
1169	!
1170	!-- Calculate df/dL
1171	f_d_ol = ( - ( z_mo / ol_u ) / ( LOG( z_mo / surf%z0(m) ) &
1172	- psi_m( z_mo / ol_u ) &
1173	+ psi_m( surf%z0(m) / ol_u ) &
1174	)**3 &
1175	+ ( z_mo / ol_l ) / ( LOG( z_mo / surf%z0(m) ) &
1176	- psi_m( z_mo / ol_l ) &
1177	+ psi_m( surf%z0(m) / ol_l ) &
1178	)**3 &
1179	) / ( ol_u - ol_l )
1180	ENDIF
1181	!
1182	!-- Calculate new L
1183	surf%ol(m) = ol_m - f / f_d_ol
1184
1185	!
1186	!-- Ensure that the bulk Richardson number and the Obukhov
1187	!-- length have the same sign and ensure convergence.
1188	IF ( surf%ol(m) * ol_m < 0.0_wp ) surf%ol(m) = ol_m * 0.5_wp
1189	!
1190	!-- If unrealistic value occurs, set L to the maximum
1191	!-- value that is allowed
1192	IF ( ABS( surf%ol(m) ) > ol_max ) THEN
1193	surf%ol(m) = ol_max
1194	EXIT
1195	ENDIF
1196	!
1197	!-- Check for convergence
1198	IF ( ABS( ( surf%ol(m) - ol_m ) / surf%ol(m) ) < 1.0E-4_wp ) THEN
1199	EXIT
1200	ELSE
1201	CYCLE
1202	ENDIF
1203	!
1204	!-- Assure that Obukhov length does not become zero
1205	IF ( ABS( surf%ol(m) ) < 1E-6_wp ) &
1206	surf%ol(m) = SIGN( 10E-6_wp, surf%ol(m) )
1207
1208	ENDDO
1209	ENDDO
1210
1211	END SUBROUTINE calc_ol
1212
1213	!
1214	!-- Calculate friction velocity u*
1215	SUBROUTINE calc_us
1216
1217	IMPLICIT NONE
1218
1219	INTEGER(iwp) :: m !< loop variable over all horizontal surf elements
1220
1221	!
1222	!-- Compute u* at horizontal surfaces at the scalars' grid points
1223	IF ( .NOT. surf_vertical ) THEN
1224	!
1225	!-- Compute u* at upward-facing surfaces
1226	IF ( .NOT. downward ) THEN
1227	!$OMP PARALLEL DO PRIVATE( z_mo )
1228	!$ACC PARALLEL LOOP PRIVATE(z_mo) &
1229	!$ACC PRESENT(surf)
1230	DO m = 1, surf%ns
1231
1232	z_mo = surf%z_mo(m)
1233	!
1234	!-- Compute u* at the scalars' grid points
1235	surf%us(m) = kappa * surf%uvw_abs(m) / &
1236	( LOG( z_mo / surf%z0(m) ) &
1237	- psi_m( z_mo / surf%ol(m) ) &
1238	+ psi_m( surf%z0(m) / surf%ol(m) ) )
1239
1240	ENDDO
1241	!
1242	!-- Compute u* at downward-facing surfaces. This case, do not consider
1243	!-- any stability.
1244	ELSE
1245	!$OMP PARALLEL DO PRIVATE( z_mo )
1246	!$ACC PARALLEL LOOP PRIVATE(z_mo) &
1247	!$ACC PRESENT(surf)
1248	DO m = 1, surf%ns
1249
1250	z_mo = surf%z_mo(m)
1251	!
1252	!-- Compute u* at the scalars' grid points
1253	surf%us(m) = kappa * surf%uvw_abs(m) / LOG( z_mo / surf%z0(m) )
1254
1255	ENDDO
1256	ENDIF
1257	!
1258	!-- Compute u* at vertical surfaces at the u/v/v grid, respectively.
1259	!-- No stability is considered in this case.
1260	ELSE
1261	!$OMP PARALLEL DO PRIVATE( z_mo )
1262	!$ACC PARALLEL LOOP PRIVATE(z_mo) &
1263	!$ACC PRESENT(surf)
1264	DO m = 1, surf%ns
1265	z_mo = surf%z_mo(m)
1266
1267	surf%us(m) = kappa * surf%uvw_abs(m) / LOG( z_mo / surf%z0(m) )
1268	ENDDO
1269	ENDIF
1270
1271	END SUBROUTINE calc_us
1272
1273	!
1274	!-- Calculate potential temperature, specific humidity, and virtual potential
1275	!-- temperature at first grid level
1276	SUBROUTINE calc_pt_q
1277
1278	IMPLICIT NONE
1279
1280	INTEGER(iwp) :: m !< loop variable over all horizontal surf elements
1281
1282	!$OMP PARALLEL DO PRIVATE( i, j, k )
1283	!$ACC PARALLEL LOOP PRIVATE(i, j, k) &
1284	!$ACC PRESENT(surf, pt)
1285	DO m = 1, surf%ns
1286
1287	i = surf%i(m)
1288	j = surf%j(m)
1289	k = surf%k(m)
1290
1291	#ifndef _OPENACC
1292	IF ( bulk_cloud_model ) THEN
1293	surf%pt1(m) = pt(k,j,i) + lv_d_cp * d_exner(k) * ql(k,j,i)
1294	surf%qv1(m) = q(k,j,i) - ql(k,j,i)
1295	ELSEIF( cloud_droplets ) THEN
1296	surf%pt1(m) = pt(k,j,i) + lv_d_cp * d_exner(k) * ql(k,j,i)
1297	surf%qv1(m) = q(k,j,i)
1298	ELSE
1299	#endif
1300	surf%pt1(m) = pt(k,j,i)
1301	#ifndef _OPENACC
1302	IF ( humidity ) THEN
1303	surf%qv1(m) = q(k,j,i)
1304	ELSE
1305	#endif
1306	surf%qv1(m) = 0.0_wp
1307	#ifndef _OPENACC
1308	ENDIF
1309	ENDIF
1310
1311	IF ( humidity ) THEN
1312	surf%vpt1(m) = pt(k,j,i) * ( 1.0_wp + 0.61_wp * q(k,j,i) )
1313	ENDIF
1314	#endif
1315
1316	ENDDO
1317
1318	END SUBROUTINE calc_pt_q
1319
1320
1321	!
1322	!-- Set potential temperature at surface grid level.
1323	!-- ( only for upward-facing surfs )
1324	SUBROUTINE calc_pt_surface
1325
1326	IMPLICIT NONE
1327
1328	INTEGER(iwp) :: k_off !< index offset between surface and atmosphere grid point (-1 for upward-, +1 for downward-facing walls)
1329	INTEGER(iwp) :: m !< loop variable over all horizontal surf elements
1330
1331	k_off = surf%koff
1332	!$OMP PARALLEL DO PRIVATE( i, j, k )
1333	!$ACC PARALLEL LOOP PRIVATE(i, j, k) &
1334	!$ACC PRESENT(surf, pt)
1335	DO m = 1, surf%ns
1336
1337	i = surf%i(m)
1338	j = surf%j(m)
1339	k = surf%k(m)
1340
1341	surf%pt_surface(m) = pt(k+k_off,j,i)
1342
1343	ENDDO
1344
1345	END SUBROUTINE calc_pt_surface
1346
1347	!
1348	!-- Set mixing ratio at surface grid level. ( Only for upward-facing surfs. )
1349	SUBROUTINE calc_q_surface
1350
1351	IMPLICIT NONE
1352
1353	INTEGER(iwp) :: k_off !< index offset between surface and atmosphere grid point (-1 for upward-, +1 for downward-facing walls)
1354	INTEGER(iwp) :: m !< loop variable over all horizontal surf elements
1355
1356	k_off = surf%koff
1357	!$OMP PARALLEL DO PRIVATE( i, j, k )
1358	DO m = 1, surf%ns
1359
1360	i = surf%i(m)
1361	j = surf%j(m)
1362	k = surf%k(m)
1363
1364	surf%q_surface(m) = q(k+k_off,j,i)
1365
1366	ENDDO
1367
1368	END SUBROUTINE calc_q_surface
1369
1370	!
1371	!-- Set virtual potential temperature at surface grid level.
1372	!-- ( only for upward-facing surfs )
1373	SUBROUTINE calc_vpt_surface
1374
1375	IMPLICIT NONE
1376
1377	INTEGER(iwp) :: k_off !< index offset between surface and atmosphere grid point (-1 for upward-, +1 for downward-facing walls)
1378	INTEGER(iwp) :: m !< loop variable over all horizontal surf elements
1379
1380	k_off = surf%koff
1381	!$OMP PARALLEL DO PRIVATE( i, j, k )
1382	DO m = 1, surf%ns
1383
1384	i = surf%i(m)
1385	j = surf%j(m)
1386	k = surf%k(m)
1387
1388	surf%vpt_surface(m) = vpt(k+k_off,j,i)
1389
1390	ENDDO
1391
1392	END SUBROUTINE calc_vpt_surface
1393
1394	!
1395	!-- Calculate the other MOST scaling parameters theta, q, (qc, qr, nc, nr)
1396	SUBROUTINE calc_scaling_parameters
1397
1398	IMPLICIT NONE
1399
1400
1401	INTEGER(iwp) :: m !< loop variable over all horizontal surf elements
1402	INTEGER(iwp) :: lsp !< running index for chemical species
1403	!
1404	!-- Compute theta* at horizontal surfaces
1405	IF ( constant_heatflux .AND. .NOT. surf_vertical ) THEN
1406	!
1407	!-- For a given heat flux in the surface layer:
1408
1409	!$OMP PARALLEL DO PRIVATE( i, j, k )
1410	!$ACC PARALLEL LOOP PRIVATE(i, j, k) &
1411	!$ACC PRESENT(surf, drho_air_zw)
1412	DO m = 1, surf%ns
1413
1414	i = surf%i(m)
1415	j = surf%j(m)
1416	k = surf%k(m)
1417
1418	surf%ts(m) = -surf%shf(m) * drho_air_zw(k-1) / &
1419	( surf%us(m) + 1E-30_wp )
1420
1421	!
1422	!-- ts must be limited, because otherwise overflow may occur in case
1423	!-- of us=0 when computing ol further below
1424	IF ( surf%ts(m) < -1.05E5_wp ) surf%ts(m) = -1.0E5_wp
1425	IF ( surf%ts(m) > 1.0E5_wp ) surf%ts(m) = 1.0E5_wp
1426
1427	ENDDO
1428
1429	ELSEIF ( .NOT. surf_vertical ) THEN
1430	!
1431	!-- For a given surface temperature:
1432	IF ( large_scale_forcing .AND. lsf_surf ) THEN
1433
1434	!$OMP PARALLEL DO PRIVATE( i, j, k )
1435	DO m = 1, surf%ns
1436	i = surf%i(m)
1437	j = surf%j(m)
1438	k = surf%k(m)
1439
1440	pt(k-1,j,i) = pt_surface
1441	ENDDO
1442	ENDIF
1443
1444	!$OMP PARALLEL DO PRIVATE( z_mo )
1445	DO m = 1, surf%ns
1446
1447	z_mo = surf%z_mo(m)
1448
1449	surf%ts(m) = kappa * ( surf%pt1(m) - surf%pt_surface(m) ) &
1450	/ ( LOG( z_mo / surf%z0h(m) ) &
1451	- psi_h( z_mo / surf%ol(m) ) &
1452	+ psi_h( surf%z0h(m) / surf%ol(m) ) )
1453
1454	ENDDO
1455
1456	ENDIF
1457	!
1458	!-- Compute theta* at vertical surfaces. This is only required in case of
1459	!-- land-surface model, in order to compute aerodynamical resistance.
1460	IF ( surf_vertical ) THEN
1461	!$OMP PARALLEL DO PRIVATE( i, j )
1462	DO m = 1, surf%ns
1463
1464	i = surf%i(m)
1465	j = surf%j(m)
1466	surf%ts(m) = -surf%shf(m) / ( surf%us(m) + 1E-30_wp )
1467	!
1468	!-- ts must be limited, because otherwise overflow may occur in case
1469	!-- of us=0 when computing ol further below
1470	IF ( surf%ts(m) < -1.05E5_wp ) surf%ts(m) = -1.0E5_wp
1471	IF ( surf%ts(m) > 1.0E5_wp ) surf%ts(m) = 1.0E5_wp
1472
1473	ENDDO
1474	ENDIF
1475
1476	!
1477	!-- If required compute q* at horizontal surfaces
1478	IF ( humidity ) THEN
1479	IF ( constant_waterflux .AND. .NOT. surf_vertical ) THEN
1480	!
1481	!-- For a given water flux in the surface layer
1482	!$OMP PARALLEL DO PRIVATE( i, j, k )
1483	DO m = 1, surf%ns
1484
1485	i = surf%i(m)
1486	j = surf%j(m)
1487	k = surf%k(m)
1488	surf%qs(m) = -surf%qsws(m) * drho_air_zw(k-1) / &
1489	( surf%us(m) + 1E-30_wp )
1490
1491	ENDDO
1492
1493	ELSEIF ( .NOT. surf_vertical ) THEN
1494	coupled_run = ( coupling_mode == 'atmosphere_to_ocean' .AND. &
1495	run_coupled )
1496
1497	IF ( large_scale_forcing .AND. lsf_surf ) THEN
1498	!$OMP PARALLEL DO PRIVATE( i, j, k )
1499	DO m = 1, surf%ns
1500
1501	i = surf%i(m)
1502	j = surf%j(m)
1503	k = surf%k(m)
1504	q(k-1,j,i) = q_surface
1505
1506	ENDDO
1507	ENDIF
1508
1509	!
1510	!-- Assume saturation for atmosphere coupled to ocean (but not
1511	!-- in case of precursor runs)
1512	IF ( coupled_run ) THEN
1513	!$OMP PARALLEL DO PRIVATE( i, j, k, e_s )
1514	DO m = 1, surf%ns
1515	i = surf%i(m)
1516	j = surf%j(m)
1517	k = surf%k(m)
1518	e_s = 6.1_wp * &
1519	EXP( 0.07_wp * ( MIN(pt(k-1,j,i),pt(k,j,i)) &
1520	- 273.15_wp ) )
1521	q(k-1,j,i) = rd_d_rv * e_s / ( surface_pressure - e_s )
1522	ENDDO
1523	ENDIF
1524
1525	IF ( bulk_cloud_model .OR. cloud_droplets ) THEN
1526	!$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
1527	DO m = 1, surf%ns
1528
1529	i = surf%i(m)
1530	j = surf%j(m)
1531	k = surf%k(m)
1532
1533	z_mo = surf%z_mo(m)
1534
1535	surf%qs(m) = kappa * ( surf%qv1(m) - surf%q_surface(m) ) &
1536	/ ( LOG( z_mo / surf%z0q(m) ) &
1537	- psi_h( z_mo / surf%ol(m) ) &
1538	+ psi_h( surf%z0q(m) / &
1539	surf%ol(m) ) )
1540	ENDDO
1541	ELSE
1542	!$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
1543	DO m = 1, surf%ns
1544
1545	i = surf%i(m)
1546	j = surf%j(m)
1547	k = surf%k(m)
1548
1549	z_mo = surf%z_mo(m)
1550
1551	surf%qs(m) = kappa * ( q(k,j,i) - q(k-1,j,i) ) &
1552	/ ( LOG( z_mo / surf%z0q(m) ) &
1553	- psi_h( z_mo / surf%ol(m) ) &
1554	+ psi_h( surf%z0q(m) / &
1555	surf%ol(m) ) )
1556	ENDDO
1557	ENDIF
1558	ENDIF
1559	!
1560	!-- Compute q* at vertical surfaces
1561	IF ( surf_vertical ) THEN
1562	!$OMP PARALLEL DO PRIVATE( i, j )
1563	DO m = 1, surf%ns
1564
1565	i = surf%i(m)
1566	j = surf%j(m)
1567	surf%qs(m) = -surf%qsws(m) / ( surf%us(m) + 1E-30_wp )
1568
1569	ENDDO
1570	ENDIF
1571	ENDIF
1572
1573	!
1574	!-- If required compute s*
1575	IF ( passive_scalar ) THEN
1576	!
1577	!-- At horizontal surfaces
1578	IF ( constant_scalarflux .AND. .NOT. surf_vertical ) THEN
1579	!
1580	!-- For a given scalar flux in the surface layer
1581	!$OMP PARALLEL DO PRIVATE( i, j )
1582	DO m = 1, surf%ns
1583	i = surf%i(m)
1584	j = surf%j(m)
1585	surf%ss(m) = -surf%ssws(m) / ( surf%us(m) + 1E-30_wp )
1586	ENDDO
1587	ENDIF
1588	!
1589	!-- At vertical surfaces
1590	IF ( surf_vertical ) THEN
1591	!$OMP PARALLEL DO PRIVATE( i, j )
1592	DO m = 1, surf%ns
1593	i = surf%i(m)
1594	j = surf%j(m)
1595	surf%ss(m) = -surf%ssws(m) / ( surf%us(m) + 1E-30_wp )
1596	ENDDO
1597	ENDIF
1598	ENDIF
1599
1600	!
1601	!-- If required compute cs* (chemical species)
1602	IF ( air_chemistry ) THEN
1603	!
1604	!-- At horizontal surfaces
1605	DO lsp = 1, nvar
1606	IF ( constant_csflux(lsp) .AND. .NOT. surf_vertical ) THEN
1607	!-- For a given chemical species' flux in the surface layer
1608	!$OMP PARALLEL DO PRIVATE( i, j )
1609	DO m = 1, surf%ns
1610	i = surf%i(m)
1611	j = surf%j(m)
1612	surf%css(lsp,m) = -surf%cssws(lsp,m) / ( surf%us(m) + 1E-30_wp )
1613	ENDDO
1614	ENDIF
1615	ENDDO
1616	!
1617	!-- At vertical surfaces
1618	IF ( surf_vertical ) THEN
1619	DO lsp = 1, nvar
1620	!$OMP PARALLEL DO PRIVATE( i, j )
1621	DO m = 1, surf%ns
1622	i = surf%i(m)
1623	j = surf%j(m)
1624	surf%css(lsp,m) = -surf%cssws(lsp,m) / ( surf%us(m) + 1E-30_wp )
1625	ENDDO
1626	ENDDO
1627	ENDIF
1628	ENDIF
1629
1630	!
1631	!-- If required compute qc* and nc*
1632	IF ( bulk_cloud_model .AND. microphysics_morrison .AND. &
1633	.NOT. surf_vertical ) THEN
1634	!$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
1635	DO m = 1, surf%ns
1636
1637	i = surf%i(m)
1638	j = surf%j(m)
1639	k = surf%k(m)
1640
1641	z_mo = surf%z_mo(m)
1642
1643	surf%qcs(m) = kappa * ( qc(k,j,i) - qc(k-1,j,i) ) &
1644	/ ( LOG( z_mo / surf%z0q(m) ) &
1645	- psi_h( z_mo / surf%ol(m) ) &
1646	+ psi_h( surf%z0q(m) / surf%ol(m) ) )
1647
1648	surf%ncs(m) = kappa * ( nc(k,j,i) - nc(k-1,j,i) ) &
1649	/ ( LOG( z_mo / surf%z0q(m) ) &
1650	- psi_h( z_mo / surf%ol(m) ) &
1651	+ psi_h( surf%z0q(m) / surf%ol(m) ) )
1652	ENDDO
1653
1654	ENDIF
1655
1656	!
1657	!-- If required compute qr* and nr*
1658	IF ( bulk_cloud_model .AND. microphysics_seifert .AND. &
1659	.NOT. surf_vertical ) THEN
1660	!$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
1661	DO m = 1, surf%ns
1662
1663	i = surf%i(m)
1664	j = surf%j(m)
1665	k = surf%k(m)
1666
1667	z_mo = surf%z_mo(m)
1668
1669	surf%qrs(m) = kappa * ( qr(k,j,i) - qr(k-1,j,i) ) &
1670	/ ( LOG( z_mo / surf%z0q(m) ) &
1671	- psi_h( z_mo / surf%ol(m) ) &
1672	+ psi_h( surf%z0q(m) / surf%ol(m) ) )
1673
1674	surf%nrs(m) = kappa * ( nr(k,j,i) - nr(k-1,j,i) ) &
1675	/ ( LOG( z_mo / surf%z0q(m) ) &
1676	- psi_h( z_mo / surf%ol(m) ) &
1677	+ psi_h( surf%z0q(m) / surf%ol(m) ) )
1678	ENDDO
1679
1680	ENDIF
1681
1682	END SUBROUTINE calc_scaling_parameters
1683
1684
1685
1686	!
1687	!-- Calculate surface fluxes usws, vsws, shf, qsws, (qcsws, qrsws, ncsws, nrsws)
1688	SUBROUTINE calc_surface_fluxes
1689
1690	IMPLICIT NONE
1691
1692	INTEGER(iwp) :: m !< loop variable over all horizontal surf elements
1693	INTEGER(iwp) :: lsp !< running index for chemical species
1694
1695	REAL(wp) :: dum !< dummy to precalculate logarithm
1696	REAL(wp) :: flag_u !< flag indicating u-grid, used for calculation of horizontal momentum fluxes at vertical surfaces
1697	REAL(wp) :: flag_v !< flag indicating v-grid, used for calculation of horizontal momentum fluxes at vertical surfaces
1698	REAL(wp), DIMENSION(:), ALLOCATABLE :: u_i !< u-component interpolated onto scalar grid point, required for momentum fluxes at vertical surfaces
1699	REAL(wp), DIMENSION(:), ALLOCATABLE :: v_i !< v-component interpolated onto scalar grid point, required for momentum fluxes at vertical surfaces
1700	REAL(wp), DIMENSION(:), ALLOCATABLE :: w_i !< w-component interpolated onto scalar grid point, required for momentum fluxes at vertical surfaces
1701
1702	!
1703	!-- Calcuate surface fluxes at horizontal walls
1704	IF ( .NOT. surf_vertical ) THEN
1705	!
1706	!-- Compute u'w' for the total model domain at upward-facing surfaces.
1707	!-- First compute the corresponding component of u* and square it.
1708	IF ( .NOT. downward ) THEN
1709	!$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
1710	!$ACC PARALLEL LOOP PRIVATE(i, j, k, z_mo) &
1711	!$ACC PRESENT(surf, u, rho_air_zw)
1712	DO m = 1, surf%ns
1713
1714	i = surf%i(m)
1715	j = surf%j(m)
1716	k = surf%k(m)
1717
1718	z_mo = surf%z_mo(m)
1719
1720	surf%usws(m) = kappa * ( u(k,j,i) - u(k-1,j,i) ) &
1721	/ ( LOG( z_mo / surf%z0(m) ) &
1722	- psi_m( z_mo / surf%ol(m) ) &
1723	+ psi_m( surf%z0(m) / surf%ol(m) ) )
1724	!
1725	!-- Please note, the computation of usws is not fully accurate. Actually
1726	!-- a further interpolation of us onto the u-grid, where usws is defined,
1727	!-- is required. However, this is not done as this would require several
1728	!-- data transfers between 2D-grid and the surf-type.
1729	!-- The impact of the missing interpolation is negligible as several
1730	!-- tests had shown.
1731	!-- Same also for ol.
1732	surf%usws(m) = -surf%usws(m) * surf%us(m) * rho_air_zw(k-1)
1733
1734	ENDDO
1735	!
1736	!-- At downward-facing surfaces
1737	ELSE
1738	!$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
1739	DO m = 1, surf%ns
1740
1741	i = surf%i(m)
1742	j = surf%j(m)
1743	k = surf%k(m)
1744
1745	z_mo = surf%z_mo(m)
1746
1747	surf%usws(m) = kappa * u(k,j,i) / LOG( z_mo / surf%z0(m) )
1748	surf%usws(m) = surf%usws(m) * surf%us(m) * rho_air_zw(k)
1749
1750	ENDDO
1751	ENDIF
1752
1753	!
1754	!-- Compute v'w' for the total model domain.
1755	!-- First compute the corresponding component of u* and square it.
1756	!-- Upward-facing surfaces
1757	IF ( .NOT. downward ) THEN
1758	!$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
1759	!$ACC PARALLEL LOOP PRIVATE(i, j, k, z_mo) &
1760	!$ACC PRESENT(surf, v, rho_air_zw)
1761	DO m = 1, surf%ns
1762	i = surf%i(m)
1763	j = surf%j(m)
1764	k = surf%k(m)
1765
1766	z_mo = surf%z_mo(m)
1767
1768	surf%vsws(m) = kappa * ( v(k,j,i) - v(k-1,j,i) ) &
1769	/ ( LOG( z_mo / surf%z0(m) ) &
1770	- psi_m( z_mo / surf%ol(m) ) &
1771	+ psi_m( surf%z0(m) / surf%ol(m) ) )
1772	!
1773	!-- Please note, the computation of vsws is not fully accurate. Actually
1774	!-- a further interpolation of us onto the v-grid, where vsws is defined,
1775	!-- is required. However, this is not done as this would require several
1776	!-- data transfers between 2D-grid and the surf-type.
1777	!-- The impact of the missing interpolation is negligible as several
1778	!-- tests had shown.
1779	!-- Same also for ol.
1780	surf%vsws(m) = -surf%vsws(m) * surf%us(m) * rho_air_zw(k-1)
1781	ENDDO
1782	!
1783	!-- Downward-facing surfaces
1784	ELSE
1785	!$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
1786	DO m = 1, surf%ns
1787	i = surf%i(m)
1788	j = surf%j(m)
1789	k = surf%k(m)
1790
1791	z_mo = surf%z_mo(m)
1792
1793	surf%vsws(m) = kappa * v(k,j,i) / LOG( z_mo / surf%z0(m) )
1794	surf%vsws(m) = surf%vsws(m) * surf%us(m) * rho_air_zw(k)
1795	ENDDO
1796	ENDIF
1797	!
1798	!-- Compute the vertical kinematic heat flux
1799	IF ( .NOT. constant_heatflux .AND. ( ( time_since_reference_point&
1800	<= skip_time_do_lsm .AND. simulated_time > 0.0_wp ) .OR. &
1801	.NOT. land_surface ) .AND. .NOT. urban_surface .AND. &
1802	.NOT. downward ) THEN
1803	!$OMP PARALLEL DO PRIVATE( i, j, k )
1804	DO m = 1, surf%ns
1805	i = surf%i(m)
1806	j = surf%j(m)
1807	k = surf%k(m)
1808	surf%shf(m) = -surf%ts(m) * surf%us(m) * rho_air_zw(k-1)
1809	ENDDO
1810	ENDIF
1811	!
1812	!-- Compute the vertical water flux
1813	IF ( .NOT. constant_waterflux .AND. humidity .AND. &
1814	( ( time_since_reference_point <= skip_time_do_lsm .AND. &
1815	simulated_time > 0.0_wp ) .OR. .NOT. land_surface ) .AND. &
1816	.NOT. urban_surface .AND. .NOT. downward ) THEN
1817	!$OMP PARALLEL DO PRIVATE( i, j, k )
1818	DO m = 1, surf%ns
1819	i = surf%i(m)
1820	j = surf%j(m)
1821	k = surf%k(m)
1822	surf%qsws(m) = -surf%qs(m) * surf%us(m) * rho_air_zw(k-1)
1823	ENDDO
1824	ENDIF
1825	!
1826	!-- Compute the vertical scalar flux
1827	IF ( .NOT. constant_scalarflux .AND. passive_scalar .AND. &
1828	.NOT. downward ) THEN
1829	!$OMP PARALLEL DO PRIVATE( i, j )
1830	DO m = 1, surf%ns
1831
1832	i = surf%i(m)
1833	j = surf%j(m)
1834	surf%ssws(m) = -surf%ss(m) * surf%us(m)
1835
1836	ENDDO
1837	ENDIF
1838	!
1839	!-- Compute the vertical chemical species' flux
1840	DO lsp = 1, nvar
1841	IF ( .NOT. constant_csflux(lsp) .AND. air_chemistry .AND. &
1842	.NOT. downward ) THEN
1843	!$OMP PARALLEL DO PRIVATE( i, j )
1844	DO m = 1, surf%ns
1845	i = surf%i(m)
1846	j = surf%j(m)
1847	surf%cssws(lsp,m) = -surf%css(lsp,m) * surf%us(m)
1848	ENDDO
1849	ENDIF
1850	ENDDO
1851
1852	!
1853	!-- Compute (turbulent) fluxes of cloud water content and cloud drop conc.
1854	IF ( bulk_cloud_model .AND. microphysics_morrison .AND. &
1855	.NOT. downward) THEN
1856	!$OMP PARALLEL DO PRIVATE( i, j )
1857	DO m = 1, surf%ns
1858
1859	i = surf%i(m)
1860	j = surf%j(m)
1861
1862	surf%qcsws(m) = -surf%qcs(m) * surf%us(m)
1863	surf%ncsws(m) = -surf%ncs(m) * surf%us(m)
1864	ENDDO
1865	ENDIF
1866	!
1867	!-- Compute (turbulent) fluxes of rain water content and rain drop conc.
1868	IF ( bulk_cloud_model .AND. microphysics_seifert .AND. &
1869	.NOT. downward) THEN
1870	!$OMP PARALLEL DO PRIVATE( i, j )
1871	DO m = 1, surf%ns
1872
1873	i = surf%i(m)
1874	j = surf%j(m)
1875
1876	surf%qrsws(m) = -surf%qrs(m) * surf%us(m)
1877	surf%nrsws(m) = -surf%nrs(m) * surf%us(m)
1878	ENDDO
1879	ENDIF
1880
1881	!
1882	!-- Bottom boundary condition for the TKE.
1883	IF ( ibc_e_b == 2 ) THEN
1884	!$OMP PARALLEL DO PRIVATE( i, j, k )
1885	DO m = 1, surf%ns
1886
1887	i = surf%i(m)
1888	j = surf%j(m)
1889	k = surf%k(m)
1890
1891	e(k,j,i) = ( surf%us(m) / 0.1_wp )**2
1892	!
1893	!-- As a test: cm = 0.4
1894	! e(k,j,i) = ( us(j,i) / 0.4_wp )**2
1895	e(k-1,j,i) = e(k,j,i)
1896
1897	ENDDO
1898	ENDIF
1899	!
1900	!-- Calcuate surface fluxes at vertical surfaces. No stability is considered.
1901	ELSE
1902	!
1903	!-- Compute usvs l={0,1} and vsus l={2,3}
1904	IF ( mom_uv ) THEN
1905	!
1906	!-- Generalize computation by introducing flags. At north- and south-
1907	!-- facing surfaces u-component is used, at east- and west-facing
1908	!-- surfaces v-component is used.
1909	flag_u = MERGE( 1.0_wp, 0.0_wp, l == 0 .OR. l == 1 )
1910	flag_v = MERGE( 1.0_wp, 0.0_wp, l == 2 .OR. l == 3 )
1911	!$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
1912	DO m = 1, surf%ns
1913	i = surf%i(m)
1914	j = surf%j(m)
1915	k = surf%k(m)
1916
1917	z_mo = surf%z_mo(m)
1918
1919	surf%mom_flux_uv(m) = kappa * &
1920	( flag_u * u(k,j,i) + flag_v * v(k,j,i) ) / &
1921	LOG( z_mo / surf%z0(m) )
1922
1923	surf%mom_flux_uv(m) = &
1924	- surf%mom_flux_uv(m) * surf%us(m)
1925	ENDDO
1926	ENDIF
1927	!
1928	!-- Compute wsus l={0,1} and wsvs l={2,3}
1929	IF ( mom_w ) THEN
1930	!$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
1931	DO m = 1, surf%ns
1932	i = surf%i(m)
1933	j = surf%j(m)
1934	k = surf%k(m)
1935
1936	z_mo = surf%z_mo(m)
1937
1938	surf%mom_flux_w(m) = kappa * w(k,j,i) / LOG( z_mo / surf%z0(m) )
1939
1940	surf%mom_flux_w(m) = &
1941	- surf%mom_flux_w(m) * surf%us(m)
1942	ENDDO
1943	ENDIF
1944	!
1945	!-- Compute momentum fluxes used for subgrid-scale TKE production at
1946	!-- vertical surfaces. In constrast to the calculated momentum fluxes at
1947	!-- vertical surfaces before, which are defined on the u/v/w-grid,
1948	!-- respectively), the TKE fluxes are defined at the scalar grid.
1949	!--
1950	IF ( mom_tke ) THEN
1951	!
1952	!-- Precalculate velocity components at scalar grid point.
1953	ALLOCATE( u_i(1:surf%ns) )
1954	ALLOCATE( v_i(1:surf%ns) )
1955	ALLOCATE( w_i(1:surf%ns) )
1956
1957	IF ( l == 0 .OR. l == 1 ) THEN
1958	!$OMP PARALLEL DO PRIVATE( i, j, k )
1959	DO m = 1, surf%ns
1960	i = surf%i(m)
1961	j = surf%j(m)
1962	k = surf%k(m)
1963
1964	u_i(m) = 0.5_wp * ( u(k,j,i) + u(k,j,i+1) )
1965	v_i(m) = 0.0_wp
1966	w_i(m) = 0.5_wp * ( w(k,j,i) + w(k-1,j,i) )
1967	ENDDO
1968	ELSE
1969	!$OMP PARALLEL DO PRIVATE( i, j, k )
1970	DO m = 1, surf%ns
1971	i = surf%i(m)
1972	j = surf%j(m)
1973	k = surf%k(m)
1974
1975	u_i(m) = 0.0_wp
1976	v_i(m) = 0.5_wp * ( v(k,j,i) + v(k,j+1,i) )
1977	w_i(m) = 0.5_wp * ( w(k,j,i) + w(k-1,j,i) )
1978	ENDDO
1979	ENDIF
1980
1981	!$OMP PARALLEL DO PRIVATE( i, j, dum, z_mo )
1982	DO m = 1, surf%ns
1983	i = surf%i(m)
1984	j = surf%j(m)
1985
1986	z_mo = surf%z_mo(m)
1987
1988	dum = kappa / LOG( z_mo / surf%z0(m) )
1989	!
1990	!-- usvs (l=0,1) and vsus (l=2,3)
1991	surf%mom_flux_tke(0,m) = dum * ( u_i(m) + v_i(m) )
1992	!
1993	!-- wsvs (l=0,1) and wsus (l=2,3)
1994	surf%mom_flux_tke(1,m) = dum * w_i(m)
1995
1996	surf%mom_flux_tke(0:1,m) = &
1997	- surf%mom_flux_tke(0:1,m) * surf%us(m)
1998	ENDDO
1999	!
2000	!-- Deallocate temporary arrays
2001	DEALLOCATE( u_i )
2002	DEALLOCATE( v_i )
2003	DEALLOCATE( w_i )
2004	ENDIF
2005	ENDIF
2006
2007	END SUBROUTINE calc_surface_fluxes
2008
2009
2010	!------------------------------------------------------------------------------!
2011	! Description:
2012	! ------------
2013	!> Calculates temperature near surface (10 cm) for indoor model or 2 m
2014	!> temperature for output
2015	!------------------------------------------------------------------------------!
2016	SUBROUTINE calc_pt_near_surface ( z_char )
2017
2018	IMPLICIT NONE
2019
2020	CHARACTER (LEN = *), INTENT(IN) :: z_char !< string identifier to identify z level
2021	INTEGER(iwp) :: i, j, k, m !< running indices
2022
2023
2024	SELECT CASE ( z_char)
2025
2026
2027	CASE ( '10cm' )
2028
2029	DO m = 1, surf%ns
2030
2031	i = surf%i(m)
2032	j = surf%j(m)
2033	k = surf%k(m)
2034
2035	surf%pt_10cm(m) = surf%pt_surface(m) + surf%ts(m) / kappa &
2036	* ( LOG( 0.1_wp / surf%z0h(m) ) &
2037	- psi_h( 0.1_wp / surf%ol(m) ) &
2038	+ psi_h( surf%z0h(m) / surf%ol(m) ) )
2039
2040	ENDDO
2041
2042
2043	CASE ( '2m' )
2044
2045	DO m = 1, surf%ns
2046
2047	i = surf%i(m)
2048	j = surf%j(m)
2049	k = surf%k(m)
2050
2051	surf%pt_2m(m) = surf%pt_surface(m) + surf%ts(m) / kappa &
2052	* ( LOG( 2.0_wp / surf%z0h(m) ) &
2053	- psi_h( 2.0_wp / surf%ol(m) ) &
2054	+ psi_h( surf%z0h(m) / surf%ol(m) ) )
2055
2056	ENDDO
2057
2058
2059	END SELECT
2060
2061	END SUBROUTINE calc_pt_near_surface
2062
2063
2064	!
2065	!-- Integrated stability function for momentum
2066	FUNCTION psi_m( zeta )
2067	!$ACC ROUTINE SEQ
2068
2069	USE kinds
2070
2071	IMPLICIT NONE
2072
2073	REAL(wp) :: psi_m !< Integrated similarity function result
2074	REAL(wp) :: zeta !< Stability parameter z/L
2075	REAL(wp) :: x !< dummy variable
2076
2077	REAL(wp), PARAMETER :: a = 1.0_wp !< constant
2078	REAL(wp), PARAMETER :: b = 0.66666666666_wp !< constant
2079	REAL(wp), PARAMETER :: c = 5.0_wp !< constant
2080	REAL(wp), PARAMETER :: d = 0.35_wp !< constant
2081	REAL(wp), PARAMETER :: c_d_d = c / d !< constant
2082	REAL(wp), PARAMETER :: bc_d_d = b * c / d !< constant
2083
2084
2085	IF ( zeta < 0.0_wp ) THEN
2086	x = SQRT( SQRT( 1.0_wp - 16.0_wp * zeta ) )
2087	psi_m = pi * 0.5_wp - 2.0_wp * ATAN( x ) + LOG( ( 1.0_wp + x )**2 &
2088	* ( 1.0_wp + x*2 ) 0.125_wp )
2089	ELSE
2090
2091	psi_m = - b * ( zeta - c_d_d ) * EXP( -d * zeta ) - a * zeta &
2092	- bc_d_d
2093	!
2094	!-- Old version for stable conditions (only valid for z/L < 0.5)
2095	!-- psi_m = - 5.0_wp * zeta
2096
2097	ENDIF
2098
2099	END FUNCTION psi_m
2100
2101
2102	!
2103	!-- Integrated stability function for heat and moisture
2104	FUNCTION psi_h( zeta )
2105	!$ACC ROUTINE SEQ
2106
2107	USE kinds
2108
2109	IMPLICIT NONE
2110
2111	REAL(wp) :: psi_h !< Integrated similarity function result
2112	REAL(wp) :: zeta !< Stability parameter z/L
2113	REAL(wp) :: x !< dummy variable
2114
2115	REAL(wp), PARAMETER :: a = 1.0_wp !< constant
2116	REAL(wp), PARAMETER :: b = 0.66666666666_wp !< constant
2117	REAL(wp), PARAMETER :: c = 5.0_wp !< constant
2118	REAL(wp), PARAMETER :: d = 0.35_wp !< constant
2119	REAL(wp), PARAMETER :: c_d_d = c / d !< constant
2120	REAL(wp), PARAMETER :: bc_d_d = b * c / d !< constant
2121
2122
2123	IF ( zeta < 0.0_wp ) THEN
2124	x = SQRT( 1.0_wp - 16.0_wp * zeta )
2125	psi_h = 2.0_wp * LOG( (1.0_wp + x ) / 2.0_wp )
2126	ELSE
2127	psi_h = - b * ( zeta - c_d_d ) * EXP( -d * zeta ) - (1.0_wp &
2128	+ 0.66666666666_wp * a * zeta )**1.5_wp - bc_d_d &
2129	+ 1.0_wp
2130	!
2131	!-- Old version for stable conditions (only valid for z/L < 0.5)
2132	!-- psi_h = - 5.0_wp * zeta
2133	ENDIF
2134
2135	END FUNCTION psi_h
2136
2137
2138	!------------------------------------------------------------------------------!
2139	! Description:
2140	! ------------
2141	!> Calculates stability function for momentum
2142	!>
2143	!> @author Hauke Wurps
2144	!------------------------------------------------------------------------------!
2145	FUNCTION phi_m( zeta )
2146	!$ACC ROUTINE SEQ
2147
2148	IMPLICIT NONE
2149
2150	REAL(wp) :: phi_m !< Value of the function
2151	REAL(wp) :: zeta !< Stability parameter z/L
2152
2153	REAL(wp), PARAMETER :: a = 16.0_wp !< constant
2154	REAL(wp), PARAMETER :: c = 5.0_wp !< constant
2155
2156	IF ( zeta < 0.0_wp ) THEN
2157	phi_m = 1.0_wp / SQRT( SQRT( 1.0_wp - a * zeta ) )
2158	ELSE
2159	phi_m = 1.0_wp + c * zeta
2160	ENDIF
2161
2162	END FUNCTION phi_m
2163
2164	END MODULE surface_layer_fluxes_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |