Home

Context Navigation

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

Last change on this file since 2234 was 2233, checked in by suehring, 8 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 78.6 KB

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