Home

Context Navigation

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

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