Home

Context Navigation

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

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