Home

Context Navigation

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

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