Home

Context Navigation

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

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