Home

Context Navigation

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

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