Home

Context Navigation

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

Last change on this file since 3456 was 3361, checked in by knoop, 6 years ago
Introduced global constant rd_d_rv=0.622
Property svn:keywords set to `Id`
File size: 86.9 KB

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