Home

Context Navigation

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

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