Home

Context Navigation

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

Last change on this file since 3272 was 3271, checked in by suehring, 6 years ago
Several critical bugfixes: initialization of surface temperature at water bodies; consider special case in initialization of vertical natural surface elements
Property svn:keywords set to `Id`
File size: 86.8 KB

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