Home

Context Navigation

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

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