Home

Context Navigation

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

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